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1.  Introduction 

It  has  been  a  known  fact  that  laminated  fiber  composites  currently 
in  use  are  relatively  weak  in  resisting  impact  loads.  Great  attention 
has  been  given  to  modeling  the  dynamic  behavior  of  composites  subjected 
to  foreign  object  impacts  and  to  the  search  for  new  forms  of  composites 
that  are  capable  of  improving  the  impact-resistant  property. 

Failure  modes  in  composites  resulting  from  the  impacts  of  a  hard 
object  and  a  soft  object  are,  in  general,  quite  different.  If  the  object 
is  relatively  rigid  and  small,  then  the  contact  time  is  short  and 
extensive  damage  is  usually  confined  to  the  neighborhood  of  the  contact 
region.  How  to  quantify  the  amount  of  damage  received  by  the  composite 
in  the  impact  zone  becomes  the  central  question  in  the  hard-object 
impact  problem. 

There  are  several  major  factors  which  could  affect  the  amount  of 
damage  in  a  laminated  composite  due  to  the  impact  of  a  hard  object. 

Among  them  are  the  mass  and  approach  velocity  of  the  object,  the  bending 
rigidity  of  the  laminate,  and  the  contact  behavior  (or  the  contact  law). 
Many  researchers  have  correlated  the  impact  velocity  with  the  damage 
for  a  given  mass.  Such  relationship  between  the  damage  and  impact 
velocity  becomes  invalid  if  the  mass  of  the  striker  or  the  bending 
property  of  the  laminate  is  changed.  The  use  of  a  single  parameter 
which  could  account  for  the  combined  effect  of  the  above  mentioned 
variables  is  highly  desirable. 

Energy  dissipation  takes  place  in  the  process  of  impact  that  results 
in  damage.  It  is  thus  reasonable  to  use  this  amount  of  energy  consumed 
in  the  impact  zone  to  measure  the  degree  of  damage  in  the  target  compos¬ 
ite  beam.  There  could  be  various  damage  modes  such  as  breakage  of 
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fibers,  cracking  in  the  matrix,  delamination,  and  plastic  deformations, 
which  could  all  contribute  to  the  energy  dissipation  in  the  impact  zone. 

It  is  conceivable  that  analytical  estimates  of  the  energies  associated 
with  these  damage  modes  are  prohibitive.  A  static  indentation  test 
which  produces  the  loading-unloading  curve  may  prove  to  be  a  simple 
means  for  determining  such  damage  energy,  since  the  energy  dissipated 
during  the  loading  and  unloading  cycle  is  simply  the  area  enclosed  by 
the  curve. 

In  this  report,  results  of  the  indentation  tests  on  glass/epoxy 
and  graphite/epoxy  laminated  composites  are  presented.  The  results 
show  that  the  loading  curve  follow  the  power  law  with  a  power  index 
1.5,  which  is  identical  to  the  classical  Hertzian  contact  law.  Sub¬ 
stantial  permanent  deformations  are  observed  even  when  loaded  at  very  low 
load  levels.  The  unloading  curves  also  follow  a  power  law. 

A  high  order  beam  finite  element  is  used  for  computing  the  dynamic 
response  of  laminated  composite  beams  subjected  to  the  impact  of  an 
elastic  sphere.  This  finite  element  includes  the  classical  elastic 
Hertzian  law  of  contact  as  well  as  the  measured  contact  law.  The  computer 
program  developed  for  this  beam  finite  element  is  listed  as  Appendix  A. 

A  simple  method  has  been  developed  for  computing  the  contact  force  and 
contact  duration.  An  estimate  of  the  contact  duration  is  needed  in  the 
finite  element  program  in  selecting  a  proper  time  increment  in  the  time 
integration  procedure.  This  method  is  found  to  be  quite  accurate  except 
for  very  thin  beams. 
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2.  Indentation  Law  for  Hard  Object  Impact  of  Composites 
2.1  Hertzian  Law  of  Contact 

When  two  solid  bodies  are  in  contact,  deformation  takes  place  in 
the  contact  zone  and  the  contact  force  results.  Once  the  contact  force 
is  obtained,  conventional  methods  for  stress  analysis  can  be  used  to  find 
the  stress  distribution  in  the  bodies.  To  determine  the  contact  force  - 
indentation  relationship  often  becomes  the  most  important  step  in  analyzing 
the  contact  problem. 

The  most  famous  elastic  contact  law  was  developed  by  Hertz  [1]  for  the 
contact  of  two  spheres  of  elastic  isotropic  materials.  The  problem  was 
solved  based  on  the  theory  of  elasticity.  A  special  case  is  that  if  the 
radius  of  one  of  the  spheres  becomes  infinite,  then  the  problem  becomes  the 
contact  of  an  elastic  sphere  and  an  elastic  half  space.  The  contact  force 
F  and  the  indentation  depth  a  were  found  to  have  the  relation 

F  =  k  (2-1) 


where 


k 

3  s 


(2-2) 


In  Eq.  (2-2),  is  the  radius  of  the  sphere,  v  is  the  Poisson's  ratio, 

E  is  the  Young's  modulus,  and  the  subscripts  1  and  2  indicate  the  two 
bodies.  Equation  (2-1)  is  usually  called  the  Hertzian  law  of  contact  for 
a  sphere  on  half  space. 

The  3/2  power  law  given  by  Eq.  (2-1)  was  found  to  be  valid  by  Willis 
[2]  for  a  rigid  sphere  pressed  on  a  transversely  isotropic  half  space. 

A  modified  contact  law  with 
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(2-3) 


was  employed  by  Sun  [3]  for  a  study  on  impact  of  laminated  composites. 

In  Eq.  (2-3),  R^,  and  E^  are  the  radius,  the  Poisson's  ratio  and  the 

Young';s  modulus  of  the  isotropic  shpere,  respectively,  and  E^  is  the 
Young's  modulus  of  the  fiber-reinforced  composite  normal  to  the  impact 
plane. 

In  applying  the  classical  Hertzian  contact  law  to  the  impact  of 
laminated  fibrous  composites  we  face  several  uncertainties.  First,  the 
half  space  assumption  is  not  valid  since  the  laminates  in  use  are  of 
finite  thickness.  Second,  the  anisotropic  and  nonhomogeneous  property  of 
laminated  compostes  may  alter  the  form  of  the  law.  Third,  the  strain  rate 
effect  which  is  not  accounted  for  by  the  Hertzian  law  may  have  significant 
effect  on  the  F-a  relation.  Except  for  the  strain  rate  effect,  the  first 
two  uncertainties  are  soTvable  by  analyzing  the  exact  contact  problem 
of  a  sphere  pressed  into  a  laminated  composite  using  three-dimensional 
elasticity.  However,  experience  tells  us  that  analytical  solutions  for 
such  contact  problems  are  extremely  difficult  to  obtain  especially  if 
permanent  deformations  are  to  be  accounted  for  during  unloadings.  Since 
unloading  paths  are  particularly  important  in  our  study,  the  experimental 
approach  is  taken  to  determine  the  law  of  contact  for  composites.  However, 
in  this  study,  the  strain  rate  effect  is  still  neglected. 


2.2  Indentation  Law  for  Laminated  Composites 
2.2.1  Theoretical  Model 

In  this  study  the  general  form  for  the  indentation  law  for  laminated 
composites  is  extended  from  the  classical  Hertzian  Law.  We  issume  that 


for  loading 
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r  «  n 

F  =  k  a 


(2-4) 


where  k  and  n  will  be  determined  experimentally.  It  is  obvious  that  when 
n  =  3/2  and  k  is  given  by  Eq.  (2-2),  this  relation  becomes  the  Hertzian 
law  for  isotropic  bodies.  It  is  noted  that  the  constant  k  has  a  very 
strange  unit  if  n  is  not  an  integer.  Also,  the  value  of  k  depends  on 
the  unit  used  for  a,  A  more  physically  meaningful  expression  may  be 
derived  by  using  a  nondimens ional  indentation  depth 

B  =  a/Rg  (2-5) 

with  which  the  indentation  can  be  written  as 

F  =  k*  3*^  (2-6) 


In  Eq.  (2-6),  k*  has  the  unit  of  force.  For  the  Hertzian  law. 


2 

s 


(2-7) 


Permanent  indentations  in  composite  targets  are  usually  generated 
even  at  relatively  low  projectile  impact  speeds.  The  permanent  indentation 
accounts  for  the  major  part  of  the  energy  loss  of  the  projectile.  Some 
energy  imparted  from  the  projectile  to  the  target  can  be  stored  in  the 
form  of  vibrational  energy  in  the  target.  As  far  as  the  local  damage 


at  the  impact  zone  is  concerned,  the  permanent  indentation  is  of  more 
interest  to  us.  For  this  reason,  the  force- indentation  law  for  the  recovery 
process  must  be  established.  In  this  study,  we  assume,  in  the  recovery 
process. 


F  =  F 


a  -  a. 


-rq 


“m"  “o 


(2-8) 


6 


where  F  is  the  maximum  contact  force  just  before  unloading  takes  place, 
m 

a  is  the  identation  corresponding  to  F  ,  and  a  is  the  permanent  indentation, 
m  mo 

This  recovery  law  was  proposed  by  Barnhart  and  Goldsmith  [4]  for  impact 
of  a  steel  ball  onto  an  armor  plate. 

2.2.2  Experimental  Results 

The  experimental  set-up  is  depicted  by  the  sketch  in  Fig.  2.1. 

The  indentation  was  measured  by  a  dial  gage  that  permits  reading  up  to 
1/5000  in.  The  dial  gage  was  mounted  on  the  loading  piston  so  that  only 
the  relative  displacement  between  the  indentor  and  the  beam  was  recorded. 

The  indentor  was  a  steel  ball  of  ^  in.  diameter.  The  beam  was  clamped  at 
both  ends  with  various  spans. 

Two  types  of  laminated  composites  have  been  tested,  namely  glass/ 
epoxy  and  graphite/epoxy.  The  glass/epoxy  was  Scotch  Ply  1002  by  the  3M 
Company.  It  contained  10  0°- plies  and  9  90°- plies  which  alternate  in  the 
layup  with  one  0°-ply  on  top  and  one  at  the  bottom.  The  thickness  of 
the  beam  was  0.19  in.  and  the  width  was  1.5  in.  The  graphite/epoxy 
specimens  were  [0/(+^45)2/02/jt^5]^  laminates.  Three  different  spans,  2  in., 

4  in.  and  6  in.,  were  used  for  the  glass/epoxy  laminates  and  two  spans, 

2  in.  and  4  in.,  were  used  for  the  graphite/epoxy  laminates. 

The  Loading  Curve 

For  the  glass/epoxy  laminate,  three  sets  of  loading  data  were  obtained 

for  each  span.  These  data  were  used  to  determine  the  best  fit  for  the 
power  law.  Eg.  (2-4),  using  the  least  squares  method.  The  results  were 

presented  in  Figs.  2. 2-2. 4.  The  power  indices  for  the  three  cases  appear  to 
be  rather  close  to  that  of  the  classical  Hertzian  law  for  isotropic  media, 
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i.e.,  n  =  1.5.  The  small  deviation  from  n  =  1.5  could  be  due  to  measure¬ 
ment  errors.  For  this  reason,  we  set  n  =  1.5  and  then  determined  k  by 
using  the  least  square  fit.  The  resulting  curves  are  shown  in  Figs. 

2. 5-2. 7.  These  curves  seem  to  fit  the  data  very  well  also. 

The  results  of  the  indentation  test  on  the  graphite/epoxy  laminated 
beams  are  presented  in  Figs.  2.8-2.10.  For  the  2-inch  span,  the  best 
least  square  fit  is  n  =  1.5;  and  for  the  4-inch  span  as  shown  in  Fig.  2.10, 
n  =  1 .5  also  yields  a  very  good  fit. 

Table  1  summarizes  the  indentation  laws  (the  loading  portion)  obtained 
from  the  experimental  results  for  a  glass/epoxy  composite  and  a  graphite/ 
epoxy  composite.  It  is  interesting  to  note  that  with  n  =  1.5,  the  values 
of  k  for  different  spans  are  almost  a  constant.  This  indicates  that  the 
indentation  law  is  ..independent  of  span.  In  other  words,  the  bending 
stress  does  not  influence  the  "contact  rigidity". 

Table  2  presents  the  indentation  laws  in  terms  of  B  and  k*  with 
n  =  1 .5  (see  Eq.  (2-6) ). 

The  Unloading  Curve 

Form  the  test  results  we  have  observed  that  permanent  deformation 

would  occur  after  an  indentation  test  no  matter  how.  small  the  load  was. 

The  unloading  paths  are  very  different  from  the  loading  path  as  can  be 

seen  from  Figs.  2.11-2.14  for  both  glass/epoxy  and  graphite/epoxy.  The 

unloading  curve  is  modeled  by  using  Eq.  (2-8)  in  which  q  and  a  have  to 

0 

be  determined.  Since  the  permanent  indentation  depth  is  difficult  to 
measure,  the  whole  data  for  each  unloading  path  were  taken  to  determine 
the  two  parameters  q  and  a^.  The  value  q  =  2.5  seems  to  yield  the  best 
overall  fit  as  shown  in  Figs.  2,11-2.14.  For  q  =  3.0  (see  Figs.  2.15-2.18) 
becomes  negative  in  some  cases. 


Table  2.  Indentation  law 
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Fig.  2.2.  Least-square  fit  of  the  contact  force  -  indentation 
relation  for  glass/epoxy  with  2-inch  span. 
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Fig.  2.3. 


Least-square  fit  of  the  contact  force-indentation 
relation  for  glass/epoxy  with  4- inch  span. 
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Fig.  2.5.  Least-square  fit  with  n  =  1.5  for  glass/epoxy  with 
2-inch  span. 
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Fig.  2.6.  Least-square  fit  with  n  =  1.5  for  glass/epoxy  with 
4- inch  span. 
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Fig.  2.7.  Least-square  fit  with  n  =  1 .5  for  glass/epoxy  with 
6- inch  span. 
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Fig.  2.8.  Least-square  fit  of  the  contact  force  -  indentation 
relation  for  graphite/epoxy  with  2-inch  span. 
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Fig.  2.9. 


Least-square  fit  of  the  contact  force  -  indentation 
relation  for  graphite/epoxy  with  4-inch  span. 
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Fig.  2,10.  Least-square  fit  with  n  =  1.5  for  graphite/epoxy 
with  4- inch  span. 
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Fig.  2.12.  Unloading  curves  for  glass/epoxy  with  4-inch  span. 
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Fig.  2.14.  Unloading  curves  for  graphite/epoxy  with  2- inch  span. 
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Fig.  2.16.  Unloading  curves  for  glass/epoxy  with  4-inch  span 


Fig.  2.17.  Unloading  curves  for  glass/epoxy 
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Fig.  2.18.  Unloading  curves  for  graphite/epoxy  with  2-inch  span. 
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3.  IMPACT  RESPONSES  BY  FINITE  ELEMENT  ANALYSIS 
3.1  The  Finite  Element 

A  beam  finite  element  with  six  degrees  of  freedom  has  been  developed 
for  the  dynamic  response  of  elastic  isotropic  beams  subjected  to  impulsive 
loadings  [5].  This  high  order  beam  element  has  been  shown  to  be  more 
efficient  than  the  conventional  element  with  four  degrees  of  freedom. 

The  element  displacement  function  is  taken  as 


2  3  4  5 

+  a2X  +  a^x  +  a^x  +  a^x  +  agX 


(3-1) 


where  v  is  the  transverse  displacement  and  a^  are  constant  coefficients. 

The  three  degrees  of  freedom  at  each  node  are  the  transverse  displacement 
V,  the  rotation  e,  and  the  curvature  k.  The  coefficients  a^  in  Eq.  (3-1) 
can  be  replaced  by  the  six  generalized  nodal  displacements  at  the  two  end 
nodes  and,  as  a  result,  the  displacement  function  can  be  alternatively 
expressed  in  terms  of  the  nodal  displacements. 

The  stiffness  and  mass  matrices  corresponding  to  the  element  displacement 
function  has  been  presented  elsewhere  [5]  and  are  reproduced  in  the  following: 
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where  Ej^I  is  the  beam  bending  rigidity,  L  is  the  length,  p  is  the  mass 
density,  and  A  is  the  cross-sectional  area.  If  the  finite  element  is 
to  be  used  for  the  analysis  of  lamianted  composite  beams,  then  the  bending 
rigidity  Ej^I  has  to  be  replaced  by  the  equivalent  bending  rigidity  D. 


(3-3) 


3.2  Impact  Response 

Based  upon  the  stiffness  and  mass  matrices  given  by  Eqs.  (3-2)  and  (3-3), 
respectively,  a  computer  program  has  been  written  specifically  for  the 
dynamic  response  of  a  beam  subjected  to  transverse  impact  of  an  elastic 
sphere.  A  finite  difference  scheme  suggested  by  Wilson  and  Clough  [6] 
was  used  to  integrate  the  time  variable  in  the  equations  of  motion.  In 
[5],  the  classical  Hertzian  law  of  elastic  contact  was  used  to  solve  a 
few  example  problems  and  excellent  results  were  found. 

The  finite  element  program  has  been  modified  for  the  analysis  of  impact 
of  laminated  beams.  The  Hertzian  indentation  laws,  Eqs.  (2-1)  with  Eq.  (2-2) 
or  Eq.  (2-3),  as  well  as  the  measured  indentation  formulas  can  be  chosen 
for  the  analysis.  Both  elastic  loading  and  actual  loading  paths  can  be 
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incorporated  in  the  program.  The  computer  program  with  a  brief  user  s 
instructions  is  presented  in  Appendix  A. 

Figures  3.1  -  3.4  show  results  for  some  example  problems  of  simply- 
supported  steel  beams,  subjected  to  impact  of  a  steel  ball.  The  diameter 
of  the  ball  is  ^  in.  The  classical  Hertzian  law  of  contact  was  used  in 
the  computation.  The  material  constants  used  are  given  by  Eq.  (4-31). 

From  these  results  it  can  be  seen  that  the  impact  velocity  has  a  great 
effect  on  the  maximum  contact  force  and  contact  duration.  The  thickness 
of  the  beam  has  little  effect  for  the  two  beam  depths  studied. 

As  reported  in  Section  2,  a  contact  of  the  steel  ball  and  the  glass/ 
epoxy  and  graphite/epoxy  composite  always  results  in  a  permanent  deformation. 
The  unloading  paths  are  substantially  different  from  the  loading  path.  If 
the  actual  unloading  paths  are  used,  the  contact  force  is  certainly  expected 
to  deviate  from  that  obtained  by  following  elastic  unloadings. 

Figures  3.5  and  3,6  present  the  results  for  a  glass/epoxy  laminated 
composite  beam  with  the  dimension  0.19  in.  D  x  1 .0  in,  W  x  7.5  in.  L. 

This  is  the  composite  beam  used  for  the  indentation  test.  The  actual 

c 

indentation  law  with  k  =  4,62  x  10  and  n  =  1.5  for  loading  and  q  =  2.5 
for  unloading  was  used  for  the  computation.  Note  that,  in  this  case,  the 
steel  ball  has  a  diameter  of  1/4  in.  same  as  that  of  the  identor  in  the 
static  indentation  test.  The  material  constants  for  composite  are 

E|_  =  5.7  X  10^  psi 
Ej  =  1 .2  X  10®  psi 
=  0.6  X  10®  psi 

^LT  ~ 

p  =  0,002016  slug/in^  (0.000168  Ib-sec^/in'^) 
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From  the  results  in  these  figures  it  can  be  seen  that  the  contact  force 
drops  more  rapidly  after  reaching  its  maximum  value  when  the  inelastic 
unloading  path  is  followed.  However,  the  total  contact  duration  does 
not  seem  to  be  affected  by  the  inelastic  unloading. 

The  finite  element  program  developed  here  can  also  be  used  in  conjunction 
with  the  experimentally  obtained  contact  law  to  compute  the  dynamic  strain 
at  any  point  on  the  beam.  The  dynamic  strain  can  be  experimentally  measured 
by  using  a  strain  gage.  By  comparing  the  measured  strain  and  that  predicted 
by  the  finite  element  solution,  it  may  be  possible  for  us  to  determine  the 
effect  of  a  result  of  this  comparison.  The  static  indentation  law  may  be 
modified  to  account  for  the  strain  rate  effect.  The  result  of  the  comparison 
will  be  reported  in  the  future. 
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Fig.  3.2.  Response  of  a  simply-supported  steel  beam 
(0.5"W  X  0.5"D  X  30"L)  subjected  to  impact 
of  a  steel  ball  with  initial  velocity  1200 
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4.  A  Simple  Method  for  Computing  Contact  Force  and  Duration  in 

Elastic  Impact 

In  using  the  finite  element  program  described  in  Section  3,  we  have 
to  choose  a  proper  time  increment  At  and  the  total  length  of  time  integration 
prior  to  the  solution.  A  poor  choice  of  At  may  result  in  poor  finite 
difference  solutions.  A  simple  way  to  obtain  an  approximate  impact  duration 
prior  to  the  use  of  the  finite  element  program  certainly  will  avoid  futile 
trials.  In  the  following,  a  simple  method  is  developed  for  computing 
an  approximate  contact  force  and  the  contact  duration. 


4.1  Impact  of  an  Elastic  Sphere  on  a  Mass  with  a  Flat  Surface 

A  simple  analysis  for  a  spherical  projectile  impacting  an  elastic 
mass  with  a  flat  surface  was  proposed  by  Timoshenko  [7]  as  follows. 

Denoting  the  mass  and  velocity  of  the  target  by  m^  and  v^,  respectively,  and 
the  mass  and  the  velocity  of  the  sphere  by  m^  and  v^,  respectively,  the 
rates  of  change  of  velocity  during  impact  are 


(4-1) 


(4-2) 


where  F  is  the  contact  force.  The  velocity  of  the  relative  approach  a 
(the  indentation)  is 


a  =  Vs  -  Vt 


(4-3) 


From  Eqs .  (4-1)  to  (4-3),  we  obtain 


a  =  -  F 


+  m^ 
t  s 


(4-4) 
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Substituting  the  Hertz  law  of  contact,  Eq.  (2-1),  in  Eq.  (4-4), 
we  obtain 


a  =  -  k^ct 


3/2 


(4-5) 


where 


m. 


(4-6) 


Integrating  Eq.  (4-5),  we  have 


1 

2  <“  -  "I* 


I 


(4-7) 


The  maximum  value  of  a,  a  ,  occurs  at  a  =  o.  We  obtain 

max 


,  5  ^  >2/5 
“max  ^4  k^  ' 


(4-8) 


This  together  with  the  Hertzian  law  yields  the  maximum  contact  force. 
From  Eq.  (4-7),  the  following  relation  is  derived; 


dt  = 


da 


,2  4  .  5/2>l/2 

(v^.  -  g  kCa  '  )  ' 


(4-9) 


By  introducing 


max 


(4-10) 


Equation  (4-9)  can  be  rewritten  as 


dt  = 


dn 


max 


(4-11) 


If  we  assume  that  the  maximum  indentation,  a  ,  is  achieved  half  way 

max 

through  the  entire  contact,  then  the  duration  of  impact  is  obtained  from 
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integrating  Eq.  (4-11)  as 


2a  fl 

max 


dn 


0 


=  2.94 


max 


(4-12) 


4.2  Equivalent  Mass  Model 

In  view  of  the  simple  formulas  given  by  Eqs.  (4-8)  and  (4-12), 
we  will  attempt  to  find  an  equivalent  mass  m^  to  represent  an  actual  beam 
or  plate.  Once  this  is  accomplished,  the  maximum  contact  force  and  the 
contact  duration  can  be  estimated  easily. 

The  equivalent  system  is  developed  based  upon  the  condition  that 
it  stores  the  same  amount  of  kinetic  and  strain  energies  as  in  the  actual 
system.  It  is  assumed  that  in  both  systems  the  strain  energies  in  the 
impactors  are  negligible  and  that  the  kinetic  energies  are  identical.  It 
is  also  assumed  that  the  spheres  do  the  same  amount  of  work  on  both  the 
actual  and  the  equivalent  targets.  With  these  assumptions,  we  conclude 
that  the  total  kinetic  energy  of  the  equivalent  mass,  K^,  should  be  equal 
to  the  kinetic  energy  K  plus  the  strain  energy  U  of  the  actual  elastic 
target,  i.e., 

K  +  U  =  (4-13) 

The  kinetic  energy  in  the  equivalent  target  system  is  simply 

From  Eq.  (4-1),  the  velocity  of  the  equivalent  mass  can  be  obtained 


by  integration  as 
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F(t)  dx 


(4-15) 


From  all  the  previous  studies,  the  contact  force  history  resembles 
a  sine  function.  In  view  of  this,  we  approximate  the  contact  force  as 
fol 1 ows 


(4-16) 


Substituting  Eq.  (3-16)  into  Eq.  (3-15)  and  integrating  from  t=o  to  t  = 
T/2  we  obtain  the  velocity  of  the  equivalent  target  at  t  =  T/2  as 


V  =  -  ^  I  F 

''t  m^  IT  ^max 


(4-17) 


Substitution  of  Eq.  (4-17)  into  Eq.  (4-14)  and  then  into  Eq.  (4-13)  lead  to 


-  1  —  (-1^ 


(4-18) 


Since  the  deflection  of  the  beam  is  proportional  to  the  applied  force  F, 

2 

both  U  and  K  contain  F  terms  and  can  be  factored  out  as 

max 


U  =  F^  U*  ,  K  =  F^  K* 
max  max 


(4-19) 


in  which  U*  and  K*  do  not  depend  on  Equation  (3-18)  can  now  be 

written  as 


2mtTr 


2  =  (U*  + 


/2 


(4-20) 


From  Eqs.  (4-6),  and  (4-8)  and  (4-12),  we  note  that  the  contact 
duration  T  is  a  function  of  the  equivalent  mass  m^.  Thus,  Eq.  (4-20)  is 
basically  a  nonlinear  equation  for  m^.  Numerical  methods  will  be  used  to 
find  solutions  for  this  equation. 
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4.3  Simply-Supported  Beam 

Consider  a  beam  of  cross-sectional  area  A  and  bending  rigidity  D. 
The  equation  of  motion  is 


3^ 

3x4 


+  pA 


3^w 


q(x,t) 


(4-21) 


where  p  is  the  average  mass  density  (over  the  thickness)  and  q(x,t)  is  a 
time  dependent  forcing  function.  For  a  homogeneous  elastic  beam,  we  have 


D  =  El 


(4-22) 


For  laminated  composite  beams,  D  is  estimated  according  to  Eq.  (4-36). 

If  the  force  is  a  concentrated  force  F(t)  applied  at  x=c,  then  the 
solution  for  Eq.  (4-21)  can  be  expressed  as  [8] 


w(x^t) 


“  W„(x)w^(c) 

I  "  " 


w^(x)dx 

0 


j-t 

F(T)sin  a)^(t-T)dT 

^0 


(4-23) 


where  w^(x)  is  the  shape  function  for  the  nth  natural  mode  of  vibration, 
and  is  the  corresponding  natural  frequency. 

For  a  simply-supported  beam,  we  obtain 


w^(x)  =  sin 


niTX 

L 


(4-24) 


and 


2 


riTT\4  D_ 
pA 


(4-25) 
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If  the  concentrated  force  is  given  by  Eq.  (4-16),  then  the  beam  deflection 
w  can  be  obtained  from  Eq.  (4-23)  as 


3 

2  F  L. 

w(x,t)  =  — 1  W  (c)  ^ 
/d  n=l  "  n^ 


4T  f  •  4. 

1-^  (sin  .ft 


4T"  -  T 


sinco^t) 


w  (x)  ,  for  t<  T 
n  - 


(4-26) 


In  Eq.  (4-26). 


T„ 

n 


(4-27) 


is  the  period  for  the  nth  mode.  The  strain  energy  and  the  kinetic  energy 
can  be  computed  in  a  straightforward  manner.  We  obtain  at  t  =  T/2 


u*  -  161-V  f 


.  16p  „ 

6^2  9l 

TT  D 


(4-28) 


where 


^1 

'  n=l 


4n^T^-T^ 


2n^T 


P  (»)l  T 

sin(n^  ) 


w„(c) 


(4-29) 


91  =  I 

'  n=l 


2  Wn  T 

-rj-j  cos(n  — )  w„(c) 


1 


4n  T  -T 


1 


(4-30) 


From  the  numerical  examples,  it  has  been  observed  that  use  of  fifty 
terms  in  the  series  in  Eqs.  (4-29)  and  (4-30)  should  provide  a  converged 
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solution.  In  all  examples  presented  in  this  section,  the  classical 
Hertzian  law  is  used. 

As  a  first  evaluation  of  the  equivalent  mass  concept,  we  consider 
a  problem  solved  by  Timoshenko  [9]  using  a  numerical  procedure  to  solve 
a  nonlinear  integral  equation.  The  steel  beam  considered  has  a  0.39  in. 

X  0.39  in.  (1  cm  X  1  cm)  cross-section  and  6.04  in.  (15.35  cm)  length. 

The  beam  is  simply-supported  at  two  ends  and  subjected  to  impact  of  a  steel 
ball  with  0.79  in.  (2  cm)  diameter.  The  material  properties  are 

E  =  31  X  10^  psi 

V  =  0.29  (4-31) 

P  =  0.00894  slug/in^  (0.000745  Ib-sec^/in^) 

It  should  be  pointed  out  that  in  the  numerical  computation,  the  value  for 
the  mass  density  as  given  by  Eq.  (4-31)  should  be  divided  by  a  factor  of 
12  if  the  length  is  given  in  inches. 

Figure  4.1  shows  the  contact  force  histories  according  to  Timoshenko's 
solution  and  the  equivalent  mass  model.  Excellent  agreement  is  noted. 

Figures  4.2  and  4.3  show  the  contact  forces  of  a  simply-supported  steel 
beam  subjected  to  impact  of  a  steel  ball  of  1/2  in.  diameter  with  different 
velocities.  The  beam  has  a  1/2  in.  x  1/2  in.  cross-section  and  is  30  in. 
long.  Both  the  equivalent  mass  model  results  and  the  finite  element  results 
are  found  to  have  a  very  close  agreement. 

The  results  for  a  thicker  steel  beam  (1/2  in.  W  x  3  in.  0  x  30  in.  L) 
with  simple  supports  are  presented  in  Figs.  4.4  and  4.5  for  v^.  =  12  in/sec.  and 
1200  in. /sec.,  respectively.  Again,  the  equivalent  mass  model  works  quite 
well  in  predicting  the  magnitude  and  duration  of  the  contact  force. 
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Figure  4.6,  shows  the  results  for  a  simple-supported  thin  steel  beam 
(0.5  in,  W  X  0.08  in.  D  x  15  in.  L)  subjected  to  the  impact  of  a  steel  ball 
of  0.5  in.  diameter  with  v.  =  100  in. /sec.  The  equivalent  mass  model  is 
able  to  predict  the  maximum  contact  force  but  not  the  contact  duration  due 
to  the  long  tail  portion. 

Figure  4.7  shows  the  contact  force  history  for  a  composite  beam  of  the  same 
dimension  and  impact  condition  as  the  previous  problem.  The  laminated 
beam  consists  of  16  piles  of  graphite/epoxy  composite.  The  ply-thickness 
is  0.005  in.  and  the  lay-up  is  (0/90/0/90)2^.  material  constants  are 

E|_  =  30  X  10^  psi,  E.^  =  0.75  x  10^  psi 

=  0.4  X  10®  psi,  =  0.25,  (4-32) 

p  =  0.00178  slug/in^  (0.000148  Ib-sec^/in^) 

The  modified  Hertzian  law  of  contact  given  by  Eq.  (2-3)  was  used  for  the 
solution.  Again,  from  Fig.  4.7  we  find  that  the  eqivalent  mass  model  is 
excellent  in  predicting  the  maximum  contact  force  but  poor  in  estimating 
the  total  contact  time.  From  the  numerical  examples  carried  out,  it  seems 
that  the  equivalent  mass  model  can  not  yield  accurate  contact  time  if  the 
target  is  too  thin. 

4.4  Simply-Supported  Rectangular  Plate 

The  plate  theory  developed  by  Whitney  and  Pagano  [10]  for  laminated 
composites  is  used  for  the  analysis.  This  plate  theory  takes  the  transverse 
shear  deformation  into  account  and  has  been  shown  by  Sun  and  Lai  [11]  to 
be  adequate  for  wave  propagation.  For  simplicity,  only  cross-ply  laminated 
plates  are  considered,  for  which  the  equations  of  motion  are  given  by 
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'’ll'''x,xx  *  “66  ''■x.yy  (“l2  +  “ee’fy.xy  '  "  ''ss'l’x  '  "  ''55''>x  =  '’^♦x 


"*12  ■"  °66>1'x,xy  “66  *y,xx  *  “22  '''y.yy  '  *44^  '  "  ''44«-„  =  P>*„  ('>-3'') 


344".y  - 


^SS’^'x.x  ^  ^55'^’xx  ^44'^y,y  ^  ^  ^44^’yy  +  P  =  phw  (4-35) 

where  a  comma  indicates  partial  differentiation,  q  is  the  lateral  load, 

w  is  the  transverse  displacement,  and  ip  are  the  rotations  of  the  plane 
2  ^ 

sectionSj  :<(=7t  /12)  is  a  shear  correction  factor,  p  is  the  average  mass 

density  (over  the  thickness),  h  is  the  piate  thickness.  I  is  the  rotary  inertia, 


1J  1J‘ 


^-n)  =  Q,,.{l,2^)dZ 


(4-36) 


In  Eq.  (4-36),  Q.^.  are  the  reduced  stiffnesses  for  the  composite  material 
For  an  isotropic  elastic  plate,  the  following  relations  exist: 


D  =  D  =  ^  ^ 

12(l-v^) 


Di2  -  vD.|^ 


^66  “  V  °11 


l  =■.  A  =  -  " 
'll  ^22  ,  2 

I  “V 


(4-37) 


A-j  2  ■“  1 


A  =  A  - 

^44  ^55  2n^ 


The  equations  of  motion  given  by  Eqs.  (4-33)  to  (4-35)  reduce  to  those  for 
the  Mindlin's  plate  theory  [12]. 


47 


If  we  separate  the  total  displacemert  into  the  bending  part,  Wj^,  and 
that  due  to  the  transverse  shear  deformation,  w^,  then  we  have 


^  =  -  ^,x 


'i'y  =  - 


(4-38) 


w  =  w.  +  w 
b  s 


In  terms  of  W|^  and  w^,  the  equation  of  motion  can  be  written  as 


^,xxx  ^  ^  ^h6\,xyy  ^  ^  ^5^,x  =  p^^,x 


(4-39) 


(D^2  +  2‘^66^'^b,xxy  ^  ^22\,yyy  ^44'^s,y  "  P’''^b,y 


(4-40) 


^  ^55^,xx  ^  %4^,yy  ^  ^ 


(4-41) 


Combining  equations  (4-39)  with  (4-40),  we  have 


^ll'^b,xxxx  ^  ^^^12  ^  ^66^'^b,xxyy  ^  *^22  '^b,yyyy  ^  ^55'^s  ,xx 


^  ''44“s,yy  “  O’K.XX  ^  "b,yy> 


(4-42) 


Equations  (4-41)  and  (4-42)  can  also  be  expressed  in  the  form 


3^  2 

LiWb  +  L2W^  =  pi  ^  V  Wj^ 

0  t 


'-2^  ^2 


(4-43) 
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where 


h  = ‘’113*2(0,2  +  2055)  2^2  +[>22  4 


9x  9y 


2  ^2 

9  9 

^2  ~  ^55  2  *  ^44  2 

^  9x  9y 


2  2 

„2  9^^  ,  9 

V  =  — 2  +  ^ 

9X  9y 


Applying  Laplace  transform  to  equations  (4-43)  yields 


*"2^s  ^  ^  ^  ^s^ 


(4-44) 


where  w^  and  w  are  the  transformed  functions  of  w.  and  w  ,  respectively, 
and  s  is  the  Laplace  transform  parameter.  Since  the  rotatory  inertia  is 
small,  it  is  neglected  in  this  study. 

Solving  Eqs.  (4-44),  we  obtain 


(Lp  -  L, )  ph  s  +  L^jLp  w  =  (-L^  +  L 


(4-45) 


We  expand  the  displacement  w  and  the  load  q  in  terms  of  the  shape  functions 


w  (x,y)  of  the  natural  modes  of  the  plate  as 
mn  ' 


«  “  n  (t)w^„  (x,y) 


>5  ■  H  q„,„  (x.J') 


(4-46) 
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For  a  simply-supported  rectangular  plate,  the  shape  function  for  the  (m,n) 
mode  is  given  by 


Sin 


mx 

a 


(4-47) 


where  a  and  b  are  the  lateral  dimensions  of  the  plate. 
Applying  Laplace  transform  to  Eq.  (4-46)  we  obtain 


"  '  n  (s)  v<^„(x,y) 
m  n 


(4-48) 


q  =  y  y  q  (s)  w  (x,y) 
i;  t;  mn  mn 
m  n 


Substitution  of  Eqs.  (4-48)  and  (4-47)  into  Eq.  (4-45)  leads  to 


S  +  0). 


(4-49) 


mn 


where 


2  1  ^mn  ^mn 


“mn  h  C  +  E 
mn  mn 


r  =  n  +  2fD  +  2D  if— +  n  fJlILi^ 

^mn  ni  ^a  ^  ^  'b  ^  ^22  ^b  ’ 


£  =  1C  A  (— +  ic  A  f— 1^ 

Sn  ^  ^  ^  ^  ^44^b  ' 


(4-50) 


The  quantity  0)^^^  is  the  angular  natural  frequency  for  the  (m,n)  mode.  If 
the  transverse  shear  deformation  is  neglected  (i.e.  the  classical  plate  theory). 


then 
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for  0  £  t  _<  T,  For  x-j  =  a/2,  =  b/2,  Eq.  (4-55)  becomes 


w(x,y,t)  = 


4p  00  00 

max  ^  ^ 


phab 


m=l  n=l 


(sinSIi  sin2fi)4 


- y  {sin  y  t  -  y -  Sin  (J3„„t) 

^  ir  \2  '  T  To) _  mn 

(- — t)  mn 


(  ^  ) 
mn 


mn 


mir  nir 

Sin  2  sin  ^ 


(4-56) 


comparing  Eq.  (4-46)  with  (4-56)  we  find 


B  (t) 
mn  ' 


4F. 


phab 


max  1 
2 


to  1 

mn 


- - j  X  sin  ^  sin  ^  x  (sin  ^  ^ 

(.  ,2  2  2  T  mn 


(-J _ ) 

mn 


(4-57) 


The  kinetic  energy  in  the  plate  at  any  time  t  is  given  by 


K  (t)  = 


ra 


(|f)'  dxdy 


(4-58) 


Substituting  Eq.  (4-57)  into  Eq.  (4-46)  and  then  into  Eq.  (4-58)  we  obtain 


H  (t)  ,  £hab  j.  j  j2^  (j) 
m  n 


(4-59) 


By  introducing  the  stiffness,  of  the  plate  system  for  the  (m,  n) 
made,  the  total  strain  energy  can  be  formally  written  as 


^  =  1  n  "mn 


(4-60) 


m  n 
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Upon  substitution  of  Eqs.  (4-59)  and  (4-60)  into  the  Lagrangian 
equation  of  motion  we  obtain 


i  phab  (4-6 

where  0  is  the  generalized  force.  From  Eq.  (4-61)  we  obtain  the  natural 
mn 

frequency  w  for  the  (m,n)  made  as 
^  mn 

2  ^  4 

“mn  phab 


mn 


(4-62) 


from  which 


_  phab  2 
mn  4  “mn 


(4-63) 


Substituting  Eqs.  (4-63)  and  (4-57)  into  Eq.  (4-60)  we  obtain 


U  (t) 


2F 

max 

phab 


I  I 

m  n 


1 


(jO  1  /  ^  > 

mn  1  -  (7;^ — j) 


mu  Hit 

- o  X  Sin  — o  sin  X 

H  \L  L  C 


mn 


x  f  *  -  if” 

mn 


(4-64) 


With  Eqs.  (4-60)  and  (4-64),  the  quantities  U*  and  K*  in  the  equivalent 


mass  mode  can  be  obtained.  We  have 
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U*  (T/2) 


phab  3 


(4-65) 


K*  (T/2)  = 


2tt^ 


phabT 


2  ^3 


(4-66) 


where 


I  I 

m  n 


1  1 


^mn  1 


(-1 _ 


/I  ■'f  '^ninT\  .4„  niTT  rnr 

V I  -  —  sin  “2 — ;  sin  sin 

mn 


mn 


(4-67) 


93  =  n 


m  n 


t2 


1  1 


0)  T 

cos  i-f-)  sin  sin  ^ 

TT  \i:  2  2  2 


^0) 

mn 


(4-68) 


Karas  [13]  considered  the  impact  of  a  steel  ball  of  2  cm  in  diameter 
on  a  simply-supported  square  steel  plate  with  a=b=20  cm  and  h  =  0,8  cm  by 
using  the  classical  plate  theory.  The  impact  velocity  of  the  ball  was  100 
cm/sec.  The  contact  force  histories  obtained  by  Karas  and  by  using  the 
equivalent  mass  model  are  shown  in  Fig,  4,8,  It  is  seen  that  the  equivalent 
mass  model  yields  a  good  estimate  of  the  maximum  contact  force  and  contact 


duration. 


BRLL  0Ifl.=l/8  IN. 


0 


1 


Simply-supported  steel  beam  (0.5"V/  x  0 
X  30"L)  subjected  to  impact  of  a  steel 
ball  at  12  in/sec. 


FORCEILB) 


Finite  Element 


I 


Fig.  4. 


-  Equivalent  Mass 


4  Simply-supported  steel  beam  (0.5"W  x  3"D 
X  30"L)  subjected  to  impact  of  a  steel 
ball  at  12  in/sec , 


FORCECLB) 


FORCE(LB) 


FORCE(LB) 


TIME  ysec 

Fig.  4.8  Contact  force  history  for  a  simply-supported 
steel  plate  (20  cm  x  20  cm  x  0.8  cm)  subjected 
to  impact  of  a  steel  ball  (2  cm  diameter)  at 
100  cm/sec. 
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5.  Conclusions 

Static  indentation  tests  have  been  performed  to  determine  the  law 
of  contact  between  a  steel  ball  and  two  laminated  composites,  namely, 
glass/epoxy  and  graphite/epoxy.  It  has  been  found  that  the  loading 
path  followed  very  well  the  power  law 

F  -  ^  1-5 

F  =  k  a 

where  F  is  the  contact  force,  k  is  a  coefficient,  and  a  is  the  indentation 
depth.  Tests  were  conducted  with  beams  clamped  at  two  ends  with  various 
spans.  The  results  indicated  that  the  indentation  law  does  not  seem  to 
depend  on  the  span  between  the  clamps.  The  experimental  results  have 
also  revealed  that  both  composites  tested  possessed  a  pronounced  inelastic 
behavior  even  at  very  low  contact  force  levels.  The  unloading  paths 
from  various  loading  points  have  been  obtained  experimentally  and  fitted 
into  a  power  law  for  the  computational  purpose. 

An  efficient  high  order  beam  finite  element  has  been  employed  together 
with  the  classical  Hertzian  contact  law  or  the  measured  contact  law  for 
analyzing  the  impact  response.  The  finite  element  program  is  capable 
of  computing  the  contact  force,  contact  duration,  and  all  the  dynamic 
responses  in  the  laminated  composite.  A  simple  method  for  estimating  the 
contact  force  and  duration  has  been  developed  and  shown  to  be  quite 
accurate  except  for  very  thin  beams. 


63 


6.  References 


H  II 

[1]  Hertz,  H.,  "Uber  die  Beruhrung  fester  elastischer  Korper" ,  Journal 
Reine  Angle  Math  (Crelle),  Vol .  92,  1881,  p.  155. 

[2]  Willis,  J.R,,  "Hertzian  Contact  of  Anisotropic  Bodies,"  Journal  of 
Mechanics  and  Physics  of  Solids,  Vol.  14,  1966,  p.  163. 

[3]  Sun,  C.Ti,  "An  Analytical  Method  for  Evaluation  of  Impact  Damage 
Energy  of  Laminated  Composites,"  ASTM  STP617,  1977,  p.  427. 

[4]  Barnhart,  K.E.,  and  Goldsmith,  W.,  "Stresses  in  Beams  during 
Transverse  Impact,"  J.  Appl .  Mech. ,  Vol.  24,  1957,  p.  440. 

[5]  Sun,  C.T.,  and  Huang,  S.N.,  "Transverse  Impact  Problems  by  Higher 
Order  Beam  Finite  Element,"  Journal  of  Computers  and  Structures, 

Vol.  5,  pp.  297-303,  1975. 

[6]  Wilson,  E.L.  and  Clough,  R.W.,  "Dynamic  Response  by  Step  by  Step 
Matrix  Analysis,"  Symposium  on  Use  of  Computers  in  Civil  Engineering, 
October  1962. 

[7]  Timoshenko,  S.P.,  Theory  of  Elasticity.McGraw-Hill ,  New  York,  1934. 

[8]  Goldsmith,  W.,  Impact,  Edward  Arnold,  London,  1960,  p.  58. 

M 

[9]  Timoshenko,  S.  "Zur  Frage  nach  der  wirkung  eines  Stosse  auf  einer 
Balken,"  Zaitschrift  fiir  Mathematik  und  Physik,  Vol.  62,  1913,  pp.  198-209. 

[10]  Whitney,  J.M.,  and  Pagano,  N.J.,  "Shear  Deformation  in  Heterogeneous 
Anisotropic  Plates,"  J.  Applied  Mechanics,  Vol.  37,  1970,  pp.  (031- 
1036. 

[11]  Sun,  C.T.,  and  Lai,  R.Y.S.,  "Exact  and  Approximate  Analysis  of  Transient 
Wave  Propagation  in  an  Anistropic  Plate,"  AIAA  Journal ,  Vol.  12,  1974, 
pp.  1415-1417. 

[12]  Mindlin,  R.D.,  "Influence  of  Rotary  Inertia  and  Shear  on  Flexural 
Vibrations  of  Isotropic,  Elastic  Plates,"  J.  Applied  Mechanics, 

Vol.  18,  1951,  pp.  31-38. 

[13]  Karas,  K.,  "Platten  Unter  Seitlichem  Stoss,"  Ingenieur-Archiv,  Vol.  10, 
1939,  pp.  237-250. 


64 


APPENDIX  A 

A  COMPUTER  PROGRAM  FOR  FINITE  ELEMENT  ANALYSIS  OF  THE  TRANSVERSE  IMPACT 
OF  A  BEAM 

The  following  is  a  description  of  the  input  data  required  to 
analyze  the  transverse  impact  of  a  beam.  The  description  is  by  card 
sections,  and  where  applicable,  the  number  of  cards  precedes  the  name. 
The  arrangement  of  the  cards  is  shown  in  Fig.  A-1. 

Heading  Card(s)  (12,  10A7) 

One  card  is  required  for  each  problem. 

Cols.  1-2  Problem  number  (NPROB) 

3-72  Arbitrary  problem  identification  (TITLE) 

2.  1-Control  Card  (915) 

Cols.  1-5  Number  of  nodal  points  (NP) 

6-10  Number  of  elements  (NE) 

11-15  Number  of  restrained  boundary  nodes  (NB) 

16-20  Number  of  output  printing  cycles  (NTM) 

21-25  Number  of  material  types  (NMAT) 

For  isotropic  materials,  this  number  is  limited  to  24 
plus  one  for  the  sphere.  However,  for  a  laminated 
composite,  this  number  can  only  be  two. 

26-30  Output  printing  frequency  in  usec  (NDIN) 

31-35  Beam  material  type  (MATP) 

0  -  if  beam  is  isotropic 
1  -  if  beam  is  a  laminated  composite 
36-50  Number  of  nodal  data  cards  (NDC) 

Explained  later. 
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41-45  Control  for  print  of  input  data  (11) 

0  -  Input  printed  at  beginning  of  first  problem  only. 

1  -  Input  printed  for  each  new  problem. 

Input  for  the  printing  scheme  outlines  the  cycle  and  frequency  at 

which  the  output  is  printed.  The  integer,  NTM,  indicates  how  many 

times  output  is  printed  after  the  sphere  makes  contact  and  the  integer, 
NDIN,  indicates  how  much  time  elapses  between  printing  of  the  output. 

In  addition,  NDIM  is  measured  in  tenths  of  a  microsecond.  As  an  example, 
if  one  wishes  to  print  output  every  5  ysec  for  10  cycles,  then  NDIM 
equals  50  (in  psec)  and  NTM=10.  Observe  that  (NDON  x  NTRM)/10  yields 
the  time  at  which  computations  stop,  in  this  case  its  50  psec. 

3.  1  -  Dimension  Card  (3F10.0) 

Cols.  1-10  Beam  thickness  (TB) 

11-20  Beam  width  (WB) 

21-30  Sphere  radius  (R) 

4.  1  -  Nodal  Impact  Card  (I5,2F10.0) 

Cols.  1-5  Impacted  node  (NQ) 

6-15  Impact  velocity  (Q2) 

16-25  Time  increment  (DT) 

5.  Element  Type  Material  Properties  Card  (s)  (I5,5F10.0) 

1  card  per  material 

Cols.  1-5  Material  number  (IMAT) 

6-15  Longitudinal  Young's  modulus  (0RT(N,1)) 

16-25  Transverse  Young's  modulus  (0RT(N,2)) 

26-35  Shear  modulus  (0RT(N,3)) 

36-45  Poisson's  ratio  (0RT(N,4)) 

46-55  Mass  density,  p  (0RT(N,5)) 
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The  last  material  card  must  contain  the  material  properties  of 
the  impacting  sphere.  If  the  sphere  and  the  beam  possess  identical 
material  properties,  then  only  one  material  card  (NMAT  =1)  is  necessary. 

6.  1  -  Identation  Law  Card  (El 0. 3 ,2F1 0. 0) 

Cols.  1-10  Loading  coefficient  k  (STF) 

11-20  Permanent  deformation  (DISPEM) 

21-30  Unloading  power  q  (QP) 

If  the  Hertzian  law  is  used  for  loading,  set  STF  =  0.0.  If  elastic 
unloading  is  followed,  then  set  DISPEM  =  0.0  and  the  input  for  QP  will 
be  ignored. 

7.  Nodal  Data  Card(s)  (2I5,  2F10.0,  1 5) 

1  card  is  required  for  each  set  of  identical  elements. 

Cols.  1-5  Beginning  node  in  the  set  (NDl) 

6-10  Final  node  in  the  set  (ND2) 

11-20  x-position  of  beginning  node  (XI) 

21-30  x-position  of  final  node  (X2) 

31-35  Element  material  type  of  set  (IMT) 

This  input  provides  information  for  the  automatic  element  generator 
in  the  program.  Given  the  above  information  for  each  set  of  identical 
elements,  the  program  computes  the  x-position  of  each  node  and  assigns 
each  element  a  material  type  and  the  Ith  and  Jth  nodes.  The  number  of 
these  cards  is  equal  to  NDC,  which  is  input  on  the  control  card. 

NOTE:  Node  1  must  begin  at  position  x  -  0. 

8.  Boundary  Conditions  Card(s)  (215) 

1  card  per  restrained  node 

Cols.  1-5  Restrained  node  (NBC) 

6-10  Boundary  condition  code  (NFIX) 

The  boundary  condition  code  is  an  integer  containir.g  three  digits. 
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Each  digit  in  the  code  is  either  1-restrained  or  0-free.  The  ones 
digit  controls  the  curvature,  the  tens  digit  controls  the  slope,  and 
the  hundreds  digit  controls  the  displacement.  As  an  example,  if  one 
node  was  clamped,  then  the  displacement  and  slope  are  zero  and  the 
curvature  is  nonzero,  or 

v  =  0  6  =  0  Kj^O  therefore  Code  110 

NOTE:  Boundary  conditions  may  be  specified  at  any  node  with  any  code. 

8.  Number  of  layers  in  laminate  (15)  (MLAYER) 

9.  Laminate  data  (15,  F5.0,  FIO.O) 

1  card  per  layer. 

Cols.  1-5  Layer  number  (L) 

6-10  Fiber  orientation  (TH) 

11-20  Layer  thickness  (TK) 

If  a  laminated  composite  beam  is  to  be  examined  under  impact,  two 
major  alterations  in  the  program  must  be  made.  This  program  provides 
for  both,  with  the  proper  indication  on  the  control  card  (MATP  =1). 

From  the  laminate  data  given,  an  equivalent  bending  rigidity  is  computed, 
or  =  El.  In  addition,  the  contact  coefficient  in  the  Hertzian 
Contact  Law  is  also  computed  differently  for  composite  beams. 

NOTE:  If  an  isotropic  beam  is  used,  skip  Cards  8  and  9. 

10.  Termination  Card 
EXAMPLE  1 

Consider  the  impact  of  a  steel  sphere  on  a  steel  cantilever  beam. 

The  dimensions  of  the  beam  are  0.5"  W  x  0.08"  0x15"  L  and  the  sphere 
has  a  diameter  of  0.5"  in.  The  initial  velocity  of  the  sphere  is  100 
in/sec.,  with  the  point  of  impact  located  at  the  mid-point.  Numerical 
solutions  are  to  be  obtained  up  to  TOO  psec  by  using  30  finite  elements. 
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The  material  constants  used  in  this  computation  are 

E  =  30  X  lO^psi,  V  =  0.25  and  p  =  0.00880  slug/in^(0. 000733  1b  -  sec^/in^). 
Note  that  the  value  of  p  in  slug  should  be  divided  by  12  if  length  is 
measured  in  inches. 

The  sample  inputs  and  outputs  for  Example  1  and  Example  2  are 
listed  following  Fig.  A-1.  The  results  for  the  contact  force,  the  dis¬ 
placement  of  the  sphere  and  the  deflection  of  the  beam  at  the  impact  point 
are  shown  in  Fig.  A-2.  The  displacement  profiles  of  the  beam  at  the 
impact  point  are  shown  in  Fig.  A-2.  The  displacement  profiles  of  the 
beam  at  various  times  are  presented  in  Fig,  A-3, 

EXAMPLE  2 

This  example  is  identical  to  the  previous  example  except  that  the 
beam  is  now  a  laminated  composite  which  consists  of  16  layers  of 
graphite/epoxy  composite.  The  ply-thickness  is  0.005”  and  the  lay-up 
is  (0/90/0/90)2^.  The  material  constants  are 

E.|.|  =  30  X  10^  psi  E22  ~  0*75  X  10^  psi  0.^2  “  0.4  x  10^  psi 

v.,2  =  0.25  ,  p  =  0.00178  slug/in^(0, 000148  lb  -  sec^/in^) 

The  corresponding  results  are  shown  in  Figs.  A-4  and  A-5. 
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'{blank  card) 


10.08Dx0.5Wxl5L  glass/epoxy  s.s.  with  Q2=100  in/sec 


(Wain  Program  &  Subroutines)  . 


PROGRAM  MAIN  (INPUT.  OUTPUT,  PLOT,  TAPES  =  INPUT, 


PFILES,  PUT,  DIM,  X=TAPE8. 


LOAD  (LGO,  RUNLIB) 
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Sample  Output  for  Example  1 


1  0.08DX0.5WX15L  ISO.  CfiHTILEUER  UITH  02=100  IM/SEC. 


NOEftL  POINTS  31 

ELEMENTS  30 

BOUNDARY  CONDITIONS  1 

OUTPUT  LIMIT  1000 

DEGREES  OF  FREEDOM  3 

MATERIALS  1 


BEAM  THICKNESS 
BEAM  UIDTH 
SPHERE  DENSITY 
SPHERE  RADIUS 


.080 

.500 

.000733 

.250 


IMPACT  NODE  IG 

IMPACT  UELOCITY  100.0 

INTEGRATION  TIME  INCREMENT!  X  E-OG  SEC)  3.500E  08 


MATERIAL  PROPERTIES 
MAT.  NO.  El 

e:2 

G12 

U12 

1  30000000.0 

30000000.0 

11500000.0 

.300 

PERMANENT  DEFORMAT I ON (IN)  ». 000000 


NDDAL  POINTS 

X 

Y 

1 

0.000 

0.000 

2 

.500 

0.000 

3 

1.000 

0.000 

4 

1.500 

0.000 

5 

2.000 

0.000 

6 

2.500 

0.000 

7 

3.000 

0.000 

8 

3.500 

0.000 

9 

4.000 

0.000 

10 

4.500 

0.000 

11 

5.000 

0.000 

12 

5.500 

0.000 

13 

G.OOO 

0.000 

14 

G.500 

0.000 

15 

7.000 

0.000 

IG 

7.500 

0.000 

17 

8.000 

0.000 

18 

8,500 

0.000 

19 

9.000 

0,000 

20 

3,500 

0.000 

21 

10.000 

0.000 

22 

10.500 

0,000 

23 

11.000 

0.000 

24 

11,500 

0.000 

25 

12.000 

0.000 

26 

12.500 

0,000 

27 

13.000 

0.000 

23 

13.500 

0.000 

23 

14.000 

0.000 

30 

14.500 

0.000 

31 

15.000 

0.000 

RHO 

.000733 


ELEMENTS 


I 

J 

1 

1 

2 

2 

2 

3 

3 

3 

4 

4 

4 

5 

5 

5 

6 

6 

G 

7 

7 

7 

8 

8 

8 

9 

3 

3 

10 

10 

10 

11 

11 

11 

12 

12 

12 

13 

13 

13 

14 

14 

14 

15 

15 

15 

16 

IG 

IG 

17 

17 

17 

18 

18 

18 

13 

19 

19 

20 

20 

20 

21 

21 

21 

22 

22 

22 

23 

23 

23 

24 

24 

24 

25 

25 

25 

26 

2B 

26 

27 

27 

27 

28 

28 

28 

29 

29 

29 

30 

30 

30 

31 

K  MAT 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 

0  0  1 


BOUriCARY  CONDI  I  IONS 
31  110 


PRINTING  SCHEME 

1.  REPORT  OUTPUT  EUERY  5.00 

2.  TERMINATE  OUTPUT  ATIDO.OO 


MSEC 

MSEC 


TYPICAL  STIFNESS  MATRIX  OF  AN  ELEMENT 


8.7F7E+04  2.1S4E+04 

2.134E+04  7.022E+03 

5.4SGE+02  2.011E+02 

-8.777E+04  -2.1S4E+04 
2.1S4E+04  3.950E+03 

-5.48SE+02  -7.314E+01 


5.486E+02  -8.777E+04  2.1S4E+04  - 

2.011E+02  -2.194E+04  3.950E+03  - 

2.743E+01  -5.48EE+02  7.314E+01 

■5.4SSE+02  8.777E+04  -2.194E+04 

7.314E+01  -2.134E+04  7.022E+03  - 

4.571E+00  5.486E+02  -2.011E+02 


TYPICAL  MASS  MATRIX  OF  AN  ELEMENT 
5.743E-0S  4.934E-07  1.S53E-08 

4.S34E-07  5.500E-03  2.28iE-09 

1.853E-03  2.281E-09  3.91St-ll 

1.587E-0G  2.39SE-07  l.i97E-08 

-2.39SE--07  -3.517E-08  -l,719E-09 
l.i97E-08  1.719E-09  8.2G3E-11 


1.587E-0G  -2.33SE-07 
2.33SE-07  -3.517E-08 
1.1S7E-0S  -1.71SE-C9 
5.743E-0S  -4.934E-07 
-4.334E-07  5.500E-03  - 

1.853E-08  -2.281E-09 


.4ESE+02 

7.314E+01 

4.571E+00 

5.48GE+02 

2.011E+02 

2,743£+01 


1.197E-08 

1.71SE-09 

8.2S3E-11 

1.85SE-08 

2.281E-09 

S.91GE-11 
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0.08DX0.5HX15L  ISO.  CANTILEUER  WITH  02=100  IN/SEC. 


TIME  ELAPSED(MSEC) 

10.500 

FORCE (LB) 

1.684E+02 

MASS  DISPLACEMENTC IN) 

9.747E-04 

MASS  UELOCITYdN/'SEC) 

7.84SE+01 

MASS  ACCEL. 

(IN/SEC2) 

-3.509E+0S 

INDENTATION(IN) 

e.l97E-04 

NODE 

DISP 

STRAIN-XX 

STRAIN-YY 

STRESS-XX 

1 

1.670E-11 

-9.832E-11 

2.9B8E-11 . 

-2.3B8E-03 

2 

-2.051E-11 

5.50BE-10 

-I.B52E-10 

l.B52E-0e 

3 

3.951E-11 

-9.258E-10 

2.777E-10 

-2.777E-02 

4 

“3.684E-11 

3.027E-10 

-9.0aiE-ll 

9.081E-03 

5 

-1.222E-10 

5.425E-09 

-1.B28E-09 

1.B28E-01 

6 

8.400E-10 

-2.G98E-08 

8.093E-09 

-8.093E-01 

7 

*2.888E~09 

8.059E-08 

-2.418E-08 

2.418E+00 

8 

6.521E-09 

-1.550E-07 

4.e50E-08 

•-4.B50E+00 

g 

-6.891E-09 

8.597E-08 

-2.579E-08 

2.579E+00 

10 

--1.677E-08 

B.471E-07 

-1.941E-07 

1.941E+01 

11 

1.002E-“07 

-2.525E-06 

7.574E-07 

-7.574E+01 

12 

~l*803E-07 

2.5B1E-0B 

-7.G83E-07 

7.B83E+01 

13 

~4.412E’-07 

1.034E-05 

-3.101E-0B 

3.101E+02 

14 

--3.972E-07 

2.215E-0S 

-B.B45E-0B 

6.645E+02 

15 

-3.358E~05 

1.188E-04 

-3.5B5E-05 

3.5G5E+03 

16 

3.578E-04 

-7.110E-04 

2.133E-04 

-2. 133E+04 

17 

~3.358E-05 

1.188E-04 

-3.5B5E-05 

3.5B5E+03 

18 

-3.972E-07 

2.215E-05 

-G.B45E-0B 

B.B45E+02 

19 

-4.412E-07 

1.034E-05 

-3.101E-0B 

3.101E+02 

20 

-1.803E--07 

2.5B1E-0S 

-7.B83E-07 

7.B83E+01 

21 

1.002E-07 

-2.525E-0B 

7.574E-07 

-7.574E+01 

22 

-1.677E-08 

G.471E-07 

-1.941E-07 

1.941E+01 

23 

“6.891E-09 

8.597E-08 

-2.579E-08 

2.579E+00 

24 

6.521E-09 

-1.550E-07 

4.B50E-08 

-4.G50E+00 

25 

-2,888E-09 

8.059E-08 

-2.418E-08 

2.418E+00 

26 

8.400E-10 

-2.G98E-08 

B.093E-09 

-8.093E-01 

27 

-1.221E-10 

5.424E-09 

-1.627E-09 

1.B27E-01 

28 

-3.697E-11 

3.0B7E-10 

-9.201E-11 

9.201E-03 

29 

3.985E-11 

-9.352E-10 

2.80BE-10 

-2.80GE-02 

30 

-'2.116E-11 

5.B98E-10 

-1.709E-10 

1.709E-02 

31 

1.306E-21 

-2.755E-10 

8.2B4E-11 

■~8. 2B4E—03 
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0.08DX0.5WX15L  ISO.  CAtfriLEUER  WITH  QS=100  IN/SEC. 


TIME  ELftPSED(MSEC) 
FORCE(LB) 

MASS  DISPLACEMENT (IN) 
MASS  UELOCITY( IN/SEC) 
MASS  ACCEL. (IN/SEC2) 
INDENTATION(IN) 


35.000 

3.223E-01 

2.332E-03 

4.839E+01 

-6.718E+03 

1.13SE-05 


mode  disp 


STRAIN-yX  STRAIN-YY 


STRE5S-XX 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 
IS 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 


-1.B39E-08 

-4.074E-08 

1.731E-07 

8.057E-08 

-1.623E-07 

-5.161E-07 

-8.232E-07 

-1.S21E-06 

-3.188E-0B 

-1.B13E-07 

2.248E-05 

-4.545E-05 

9.28BE-05 

-3.946E-04 

1.079E-03 

2.322E-03 

1.079E-03 

-3.S4SE-04 

9.286E-05 

-4.545E-05 

2.248E-05 

-l,612E-07 

-3.188E-0B 

-1.B21E-06 

-8.225E-07 

-5.171E-07 

-1.620E-07 

8.2B3E-08 

1.B84E-07 

-3.759E-08 

-2.132E-18 


1.939E-07 
1.315E-06 
-3. 1B3E-0B 
-1.543E-06 
1.457E-06 
6.35BE-0B 
8.349E-0B 
1.321E-05 
1.493E-05 
-1.602E-05 
-5.990E-05 
1.443E-04 
-2.545E-04 
4.B3BE-04 
-B.524E-05 
-4.584E-04 
-B.524E-05 
4.B36E-04 
-2.545E-04 
1.443E-04 
-5.990E-05 
-1.602E-05 
1.493E-05 
1.322E-05 
8.333E-0B 
6.379E-06 
1.454E-06 
-1.598E-0B 
-3.0S7E-0B 
1.254E-0B 
5.192E-08 


-5.817E-08 

-3.944E-07 

9.490E-07 

4.B30E-07 

-4.371E-07 

-1.907E-0B 

-2.505E-0B 

-3.9B4E-0B 

-4.478E-06 

4.806E-06 

1.797E-05 

-4.329E-05 

7.636E-05 

-1.391E-04 

1.957E-05 

1.375E-04 

1.957E-05 

-1.391E-04 

7.B3BE-05 

-4.329E-05 

1.797E-05 

4.807E-06 

-4.478E-06 

-3.9B6E-0B 

-2.500E-06 

-1.914E-06 

-4.362E-07 

4.793E-07 

9.171E-07 

-3.7B2E-07 

-1.558E-08 


5.817E+00 

3.944E+01 

-9.490E+01 

-4.630E+01 

4.371E+01 

i.907E+02 

2.505E+02 

3.964E+02 

4.478E+02 

-4.80BE+02 

-1.797E+03 

4.329E+03 

-7.63BE+03 

1.391E+04 

-1.357E+03 

-1.37SE+04 

-1.957E+03 

1.391E+04 

-7,63BE+03 

4.329E+03 

-1.797E+03 

-4.807E+02 

4.478E+02 

3.9B6E+02 

2.500E+02 

1.914E+02 

4.3B2E+01 

-4.793E+01 

-9.171E+01 

3.7B2E+01 

1.558E+00 
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Sample  Input  for  Example  2 
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Sample  Output  for  Example  2 


1  0.08DX0.5WX15L  COMP.  CANTILEUER  WITH  02=100  IN/SEC 


NODAL  POINTS  31 

ELEMENTS  30 

BOUNDARY  CONDITIONS  1 

OUTPUT  LIMIT  1000 

DEGREES  OF  FREEDOM  3 

MATERIALS  2 


BEAM  THICKNESS 

.080 

BEAM  WIDTH 

.500 

SPHERE  DENSITY 

.000733 

SPHERE  RADIUS; 

.250 

IMPACT  NODE 

IB 

IMPACT  UELOCITY 

100.0 

l.OOOE-07 

INTEGRATION  TIME  INCREMENTC  X  E-OB  SEC) 

MATERIAL  PROPERTIES 

MAT.  NO.  El 

E2 

G12 

1  30000000.0 

750000.0 

400000.0 

2  30000000.0 

30000000.0 

11500000.0 

PERMANENT  DEFORMATIONC IN)  *.000000 


U12 

.250 

.300 


ADD  MATRIX  „  ^ 

1.232E+0B  1.502E+04  -5.912E-3S 

1.502E+04  1.232E+0S  -1.9S4E-09 

-5.91EE-39  -1.394E-09  3.200E+04 

-2.910E-11  3. 183E-12  -9.788E-55 

3.183E-12  3.402E-10  -4.653E-25 

-9.788E-55  -4.B53E-25  4.547E-12 

NODAL  POINTS 


X 

Y 

1 

0.000 

0.000 

2 

.500 

0.000 

3 

1.000 

0.000 

4 

1.500 

0.000 

5 

2.000 

0.000 

B 

2.500 

0.000 

7 

3.000 

0.000 

8 

3.500 

0.000 

8 

4.000 

0.000 

10 

4.500 

0.000 

11 

5.000 

0.000 

12 

5.500 

0.000 

13 

B.OOO 

0.000 

14 

G.500 

0.000 

15 

7.000 

0.000 

16 

7.500 

0.000 

17  ■ 

8.000 

0.000 

18 

8.500 

0.000 

19 

9.000 

0.000 

20 

9.500 

0.000 

21 

10.000 

0.000 

22 

10.500 

0.000 

23 

11.000 

0.000 

24 

11.500 

0.000 

25 

12.000 

0.000 

26 

12.500 

0.000 

27 

13.000 

0.000 

28 

13.500 

0.000 

29 

14.000 

0.000 

30 

14.500 

0.000 

31 

15.000 

0.000 

-2.910E-11  3.183E-12  -9.788E-55 

3.183E-12  3.402E-10  -4.B53E-25 

-9.788E-55  -4.653E-25  4.547E-12 

7.742E+02  8.013E+00  -2.5B2E-42 

8.013E+00  5.338E+02  -8.S39E-13 

-2.5G2E-42  -8.633E-13  1.707E+01 


RHO 

000148 

000733 


ELEMENTS 


I 

J 

K 

MAT 

1 

1 

2 

0 

0 

1 

a 

2 

3 

0 

0 

1 

3 

3 

4 

0 

0 

1 

4 

4 

5 

0 

0 

1 

5 

5 

6 

0 

0 

1 

6 

6 

7 

0 

0 

1 

7 

7 

8 

0 

0 

1 

8 

8 

9 

0 

0 

1 

9 

9 

10 

0 

0 

1 

10 

10 

11 

0 

0 

1 

11 

11 

12 

0 

0 

1 

12 

12 

13 

0 

0 

1 

13 

13 

14 

0 

0 

1 

14 

14 

15 

0 

0 

1 

15 

15 

16 

0 

0 

1 

16 

16 

17 

0 

0 

1 

17 

17 

18 

0 

0 

1 

18 

18 

19 

0 

0 

1 

19 

19 

20 

0 

0 

1 

20 

20 

21 

0 

0 

1 

21 

21 

22 

0 

0 

1 

22 

22 

23 

0 

0 

1 

23 

23 

24 

0 

0 

1 

24 

24 

25 

0 

0 

1 

25 

25 

26 

0 

0 

1 

26 

26 

27 

0 

0 

1 

27 

27 

28 

0 

0 

1 

28 

28 

29 

0 

0 

1 

29 

29 

30 

0 

0 

1 

30 

30 

31 

0 

0 

1 

BOUNDARY  CONDITIONS 
31  110 

PRINTING  SCHEME 

1.  REPORT  OUTPUT  EUERY  5.00  MSEC 
a.  TERMINATE  OUTPUT  ATIOO.OO  MSEC 

TYPICAL  STIFNESS  MATRIX  OF  AN  ELEMENT 
l.OGaE+05  a.654E+04  B.63GE+0a  -l.OGEE+OS  a.G54E+04  -G.GSGE+Oa 
a.G54E+04  8.494E+03  a.433E+0a  -a.G54E+04  4.778E+03  -8.848E+01 


6.636E+02 

2.433E+02 

3.318E+01 

-B.B3BE+0a  8.848E+01 

5.530E+00 

“1.0G2E+05 

-2.G54E+04  ■ 

-G.63GE+02 

l.OBaE+05  -a.B54E+04 

6,636E+02 

2.654E+04 

4.778E+03 

8.848E+01 

-a.B54E+04  8.494E+03 

-2.433E+02 

-G.63BE+02 

-8.848E+01 

5.530E+00 

B.G3BE+0a  -2.433E+02 

3.318E+01 

TYPICAL  MASS  MATRIX  OF 

AN  ELEMENT 

l.lGOE-06 

9.963E-08 

3.751E-09 

3.203E-07  -4.837E-08 

2,416E-09 

9.9G3E-08 

l.lllE-08 

4.G05E-10 

4.837E-08  -7.101E-09 

3.470E*-10 

3.751E-09 

4.G05E-10 

2.002E-11 

2.41BE-09  -3.470E-10 

1,GB8E-11 

3.203E-07 

4.837E-08 

2.41GE-09 

I.IGOE-OB  -9.9B3E-08 

3.751E-09 

•“4.837E-08 

-7.101E-09  ■ 

“3.470E-10 

-9.9G3E-08  l.lllE-08 

-’4.605E-10 

2.416E-09 

3,470E-10 

1.6G8E-11 

3.751E-09  -4.B05E-10 

2.002E-11 

0.08DX0.5MX15L  COMP.  CANTILEUER  WITH  02=100  IM/SEC 


TIME  ELAPSED(MSEC) 

20.000 

FORCE < LB) 

2.989E+01 

MASS 

DISPLACEMENT(IN) 

1.3B3E-03 

MASS  UELOCITY(IN/SEC) 

9.407E+01 

MASS 

ACCEL 

.  ( IN/SEC2) _ 

-6.230E+05 

IMDENTATIONCIN) 

1.557E-03 

NODE 

DISP 

STRAIM-XX 

STRAIN-YY 

STRESS-XX 

1 

4.354E-0S 

-1.203E-08 

3.008E-09 

-3.B10E-01 

2 

1.211E-09 

-1.140E-08 

2.850E-09 

-3.420E-01 

3 

1.594E-09 

-1.941E-08 

4,853E-09 

-5.824E-01 

4 

4.0B3E-09 

-5.353E-08 

1.338E-08 

-1.60BE+00 

5 

1.002E-08 

-1.048E-07 

2.619E-08' 

-3.143E+00 

6 

2.165E-08 

-1.112E-07 

2.780E-08 

-3.33BE+00 

7 

2.325E-08 

-4.028E-08 

l,007E-08 

-1.209E+00 

8 

-7.712E-08 

3.099E-07 

“7.748E-08 

9.297E+00 

9 

-2.B02E-0a 

-1.B33E-07 

4.081E-08 

-4.898E+00 

10 

3.6BBE-07 

-2.982E-07 

7.454E-08 

-8.945E+00 

11 

-1.191E-0B 

8.802E-07 

-2.201E-07 

2.B4IE+01 

12 

2.941E-0B 

-8.330E-08 

2.098E-08 

-2.517E+00 

13 

-1.8B2E-06 

-9.902E-06 

2.476E-06 

-2.971E+02 

14 

-3.526E-05 

3.434E-05 

-8.586E-06 

1.030E+03 

15 

1.289E-04 

3.357E-05 

-9.892E--06 

1.187E+03 

16 

4.115E-04 

-1.934E-04 

4,835E-05 

-5.802E+03 

17 

1.283E-04 

3.957E-05 

-9.892E“06 

1 . 187E+03 

18 

-3.52BE-05 

3.434E-05 

-8.586E-06 

1.030E+03 

15 

-1.862E-0B 

-9.902E-0B 

2.476E-06 

-2.971E+02 

20 

2.941E-0B 

-8.370E-08 

2.092E-08 

-2.511E+00 

21 

-1.191E-06 

8.798E-07 

-2.199E-07 

2.639E-I-01 

22 

3.BB6E-07 

-2.978E-07 

7.44GE-08 

-8.935E+00 

23 

-2.B03E-08 

-1.B23E-07 

4.073E-08 

-4.888E+00 

24 

-7.707E-08 

3.088E-07 

-7.721E-08 

9.2B5E+00 

25 

2.321E-08 

-3.940E-08 

9.850E-09 

-1.182E+00 

26 

2.1B4E-08 

-1.112E-07 

2.779E-08 

-3.335E+00 

27 

1.003E-08 

-1.050E-07 

2.624E-08 

-3.149E+00 

28  • 

4.071E-09 

-5.328E-08 

1.332E-08 

-1.599E+00 

29 

1.B9BE-09 

-2.162E-08 

5.406E-09 

-B.487E-01 

30 

1.13SE-09 

-9.681E-09 

2.420E-09 

-2.904E-01 

31 

2.507E-18 

-3.823E-08 

9.556E-09 

-1.147E+00 
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0.08DX0.5WX15L  COMP.  CANTILEUER  WITH  02=100  IN/SEC 


TIME  ELAPSED (MSEC) 

90.000 

FORCE (LB) 

2.785E+00 

MASS  DISPLACEMENT(IN) 

7.288E-03 

MASS  UELOCITY( IN/SEC) 

B.770E+01 

MASS  ACCEL. 

CIN/SEC3) 

-5.805E+04 

INDENTATIONdN) 

3.2B4E-04 

NODE 

DISP 

STRAIN-XX 

STRAIN-YY 

1 

3.708E-05 

-2.198E-08 

5.49GE-09 

a 

-3.087E-05 

2.453E-05 

-B. 132E-0B 

3 

2.007E-05 

-1.937E-0B 

4.844E-07, 

4 

5.413E-05 

-3.431E-05 

8.577E-0G 

5 

-7.G5BE-0S 

1.S02E-05 

-4.004E-0G 

6 

-1.137E-04 

5.4G7E-05 

-1.3G7E-05 

7 

1.3GSE-04 

-3.582E-0B 

8.954E-07 

8 

3.434E-04 

-9.G9GE-05 

2.424E-05 

9 

5.487E-0B 

-7.98GE-05 

1.99GE-05 

10 

-7.572E-04 

B.927E-05 

-1.732E-05 

11 

-l.lOOE-03 

1.993E-04 

-4.982E-05 

IS 

-2.GB4E-04 

1.998E-04 

-4.995E-05 

13 

1.748E-03 

7.519E-05 

-1.880E-05 

14 

4.217E-03 

-8.072E-05 

2.018E-05 

15 

G.202E-03 

-1.999E-04 

4.998E-05 

16 

G.9G9E-03 

-2.B14E-04 

B.53GE-05 

17 

B.202E-03 

-1.999E-04 

4.998E-05 

18 

4.217E-03 

-8.072E-05 

2.018E-05 

19 

1.748E-03 

7.519E-05 

-1.880E-05 

SO 

-2.GG4E-04 

1.998E-04 

-4.995E-05 

SI 

-l.lOOE-03 

1.993E-04 

-4.982E-05 

SS 

-7.572E-04 

G.927E-05 

-1.732E-05 

S3 

5.487E-0B 

-7.98GE-05 

1.997E-05 

S4 

3.434E-04 

'■9.B9BE-05 

2.424E-05 

S5 

1.3B5E-04 

-3.579E-0B 

8.949E-07 

S6 

-1.137E-04 

5.4G7E-05 

-1.3G7E-05 

S7 

-7.G54E-05  . 

1.603E-05 

-4.007E-0G 

S8 

5.425E-05 

-3.428E-05 

8.571E-0G 

S9 

2.033E-05 

-2.15BE-0B 

5.389E-07 

30 

-3.311E-05 

2.110E-05 

-5.275E-0G 

31 

-3.5SEE-15 

-3.097E-05 

7.743E-0S 

STRESS-XX 

-G.595E-01 

7.358E+02 

-5.812E+01 

-1.029E+03 

4.805E+02 

1.B40E+03 

-1.074E+02 

-2.909E+03 

-2.39GE+03 

2.078E+03 

5.978E+03 

5.994E+03 

e.25GE+03 

-2.422E+03 

-5.998E+03 

-7.843E+03 

-5.998E+03 

-2.422E+03 

2.25GE-)-03 

5.994E+03 

5.978E+03 

2.078E+03 

-2.39GE+03 

-2.909E+03 

-1.074E+02 

1.G40E+03 

4.808E+02 

-1.023E+03 

-B.4G7E+01 

G.330E+02 

-9.291E+02 
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Fig.  A-5  Displacement  profiles 
composite  beam. 


at  various  times  after  impact  of  the 
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Listing  of  Program 


PROGRAM  MAIN  ( INPUT, OUTPUT, PLOT, TAPE5=INPUT, TAPEe=OUTPUT, TAPEl 1 , TA 
1PE8) 


CONTROL  MAIN  PROGRAM 


COMMON  /CONTR/  TITLEC?) , NP, NE, NB, NDF, NCN, NLD, NMAT, NSZF, LI, NT4, NDIN 
1,MATP,NPR0B 

COMMON  /TIME/  T, DT, DDT, TAU, KCON, KCNT 
COMMON  /DISP/  01, Q2, 03, 010,020,030 
COMMON  /DIME/  TB,WB,PB,NQ,D11 
COMMON  /SPHERE/  STF,R,CABU(10),QKONST(10) 

COMMON  /PLASTIC/  DISPEM,NDISPEM,FORSPM,DISPM,OP 

COMMON  CORD(100,2),NOP(200,4), IMAT(200),ORT(25,5),NBC(25),NFIX(25) 
1,R1(200),R2(200),R3(200),R10(200),R20(200),R30(200),FORS(200),SM(2 
200,15),SK(200, 15), ISP(200,15),SMPEM(200,15),ESTIF(12, 1E),EMASS<12, 
312),NFIXK(25) 

COMMON  /COMP/  0BR(3,3,25),ABD(6,S),TH(25),ZK(25),MLAYER 
COMMON  /PLOT/  NN,TTC25),FF(25),W(25),U(25) 


INITIALIZE  TAPE  NO. 

AND  NUMBER  OF  CORNER  NODE  MAX. 

NT4=11 

NCN=2 

NN=1 


PROBLEM  IDENTIFICATION 
CALL  PLOTS 

101  READ  (5.108)  NPROB, (TTTLECI) . I=l»7) 

IF  (NPROB.EQ.O)  GO  TO  105 

DO  102  KG=1.200 
R10(KG)=0. 

R20CKG)=0. 

R30(KG)=0. 

102  R3(KG)=0o 

READ  INPUT  GEOMETRV  AND  PROPERTIES 

CALL  GDATA 

NDISPEM=0 

T=0. 

TAU=2. 

KCON=0 

DDT=DT*DT 

LOOP  ON  NO  OF  PROBLEMS 

REWIND  NT4 
NSZF=NP*NDF 
CALL  FORNK 
CALL  FORMM 
DO  103  LI=1.NLD 
KCNT^l 

READ  LOADS 

CALL  LOAD 

FORM  THEN  SOLUE  SIMULTANEOUS  EQUATIONS 

CALL  HMTQ 
CALL  SOLUE 
CALL  INTEGTN 

ITERATION  2 

KCNT=2 
CALL  LOAD 
CALL  HMTQ 
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A 

A 

A 

A 

A 

A 

A 

A 
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A 

A 

A 
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A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 
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A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 
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A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 

A 


2 

3 

4 

5 
G 

7 

8 
9 

10 

11 

12 

13 

14 

15 
IG 

17 

18 

19 

20 
21 
22 

23 

24 

25 
2G 

27 

28 

29 

30 

31 

32 

33 

34 

35 
3G 

37 

38 

39 

40 

41 

42 

43 

44 

45 
4G 

47 

48 

49 

50 

51 

52 

53 

54 

55 
5G 

57 

58 

59 

60 
61 
G2 

63 

64 

65 

66 

67 

68 

69 

70 

71 
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CALL  SOLUE  A  72 

CALL  INTEGTN  A  73 

T=T+DT  A  74 

IF  (T.GT.IOO.E-G)  GO  JO  104  A  75 

IF  ( LI. EQ. 10000)  GO  TO  104  .  A  7G 

103  CONTINUE  A  77 

104  WRITE  (G»10G)  A  78 

WRITE  (G,107)  ((TT(I),FF(I).W(I),U(I)),I=1,NN)  A  73 

WRITE  (8.107)  ((TT(I),FF(I),W(I).U(I)),I=1»NN)  A  80 

CALL  FACTOR  (0.8)  A  81 

CALL  PLOT  (0.0. 2. 0.3)  A  82 

CALL  SCALE  (TT.G.0.21. 1)  A  83 

CALL  SCALE  (FF.S.0.21. 1)  A  84 

CALL  SCALES  (S.O.W.Sl. 1.U.21. 1)  A  85 

A  8G 

W(22)=U(22)=TT(22)=FF(22)=0.0  A  87 

TT(23)=20.  A  88 

FF(23)=20.  A  83 

W(23)=U(23)=0.001  A  30 

A  31 

CALL  AXIS  (0.0. 0.0. 10HTIME(  SEC). -10. G.O. 0.0. TT(22).TT(23).0)  A  32 

CALL  AXIS  (0.0,0.0.3HFORCE(LB).3.3.0.30.0.FF(22).FF(23).-1)  A  33 

CALL  AXIS  (G.0.0.0.14HDISP(0.001  IN). 14.9. 0.30.0. W(22).W(23).-1)  A  34 

CALL  LINE  (TT.FF.21. 1.1.3)  A  95 

CALL  DSHLINE(TT.W.21. 0.1. 0.1.1)  A  9G 

CALL  DSHLINE  (TT.U.21. 0.05, 0.05. 1)  A  97 

CALL  PLOT  (6. 0,9. 0.3)  A  98 

CALL  PLOT  (0.0, 9. 0,2)  A  99 

CALL  SYMBOL  ( 1 . 0,9.3, 0. 1. TITLE, 0. 0.70)  A  100 

CALL  SYMBOL  (3.5. 8.5. 0. 1. 17HBALL  DIA.=l/2  IN., 0.0, 17)  A  101 

GO  TO  101  A  102 

105  CALL  PLOT  (0,0,999)  A  103 

STOP  A  104 

A  105 

1 OB  FORMAT  (1H1.4X. 10HTIME(MSEC),GX.9HF0RCE(LB),2X, 13HBALL  DISP(IN),2X  A  lOG 

1,13HBEAM  DISP(IN))  A  107 

107  FORMAT  (4E15.3)  A  108 

108  FORMAT  (I2,7A10)  A  109 

A  110 

END  A  111 

SUBROUTINE  GDATA  B  2 

COMMON  /CONTR/  TITLE(7),NP,NE,NB.NDF.NCN.NLD,NMAT.NSZF,LI.NT4.NDIN  B  3 

l.MATP.NPROB  B  4 

COMMON  /TIME/  T.DT. DDT. TAU.KCON.KCNT  B  5 

COMMON  /DISP/  Q1.Q2,Q3.Q10.Q20.Q30  B  G 

COMMON  /DIME/  TB.WB.PB.NQ.Dll  B  7 

COMMON  /SPHERE/  STF,R,CABU(10),QKONST(10)  B  8 

COMMON  /PLASTIC/  DISPEM.NDISPEM.FORSPM.DISPM.QP  B  9 

COMMON  CORD(100,2).N0P(200,4),IMAT(200),ORT(25,5),NBC(25).NFIX(25)  B  10 

l.Rl(200),R2(200).R3(200),R10(200).R20(200).R30(200).FORS(200).SM(2  B  11 

200,15),SK(200, 15),ISP(200,15),SMPEM(200.15),ESTIF(12, 12),EMASS(12,  B  12 

312),NFIXK(25)  B  13 

COMMON  /COMP/  QBR(3,3,25),ABD(G.G).TH(25).ZK(25),MLAYER  B  14 

COMMON  /PLOT/  NN.TT(25).FF(25),W(25),U(25)  B  15 

B  IG 

READ  AND  PRINT  TITLE  AND  CONTROL  B  17 

B  18 

WRITE  (GfllG)  NPROB. (TITLE(I),I=1»7)  B  19 

WRITE  (8, IIG)  NPROB, (TITLE(I), 1=1,7)  B  20 

READ  (5,10G)  NP,NE,NB,NTM,NMAT,NDIN,MATP,NDC,Ii  B  21 

NDF=3  B  22 

NLD=NDIN*NTM  B  23 

FLD=FLOAT(NLD)/10.  B  24 

FDIN=FLOAT(NDIN)/10.  B  25 

READ  (5,113)  TB,WB,R,NQ,Q2,DT  B  2G 

WRITE  (E,107)  NP,NE,NB,NLD,NDF,NNAT  B  27 

NLD=NLD+1  B  28 

B  29 

READ  AND  PRINT  MATERIAL  DATA  B  30 

SPHERE  DATAs  L=NMAT  (LAST  MAT.  CARD)  B  31 


UUCJ  U  U  U  U  tJ  U  U  L)  UUCJCJ  UCJU 


86 


C  B  32 

READ  C5.112)  (L,(0RT(L,I),I=l,5),N=l,mAT)  B  33 

PB=0RTCNMAT,5)  B  34 

WRITE  (6,122)  TB,WB,PB,R,NQ,Q2,DT  B  35 

NQ=(NQ-“1)^*3+1  B  36 

WRITE  (6,121)  B  37 

WRITE  (6,115)  (N,(DRT(N,I),I=1,5),N=1,NMAT)  B  38 

E  39 

READ  IMDENTATION  DATA  B  40 

B  41 

READ  (5*111)  STF.DISPEM.QP  B  42 

WRITE  (B*123)  DISPEM  B  43 

IF  (DISPEIi.NE.0.0)  WRITE  (G.  124)  QP  B  44 

B  45 

READ  NODAL  POINT  DATA  B  4B 

AND  B  47 

READ  ELEMENT  DATA  B  48 

B  49 

DO  102  1=1. NDC  B  50 

READ  (5,114)  ND1*ND2,X1,X2*IMT  B  51 

EL=(X2-X1)/FL0AT(ND2-ND1)  B  52 

C0RD(ND1, 1)=X1  B  53 

C0RD(ND2,1)=X2  B  54 

CORD(ND2,2)=0,()  B  55 

C0RD(ND1,2)=C0RD(ND2,2)  B  5G 

NDD=ND2-1  B  57 

DO  101  J=ND1,NDD  B  58 

CORD(J-i-l*l)=CORD(J,l)+EL  B  59 

CORD(J+1,2)=0.0  B  BO 

NOP(J,l)=J  B  B1 

N0P(J*2)=J+1  B  G2 

NOP(J,4)=0  B  G3 

N0P(J,3)=N0P(J,4)  B  64 

IMAT(J)=IMT  B  B5 

101  CONTINUE  B  GB 

102  CONTINUE  B  G7 

B  68 

READ  BOUNDARY  DATA  B  G9 

B  70 

READ  (5,110)  (NBC(I),NFIX(I),I=1,NB)  B  71 

IF  (MATP.EQ.i)  CALL  CMPD  B  72 

B  73 

ISOTROPIC  MATP=0.0  B  74 

COMPOSITE  MATP=1.0  B  75 

B  7G 

IF  (Il.NE.O)  GO  TO  103  B  77 

B  78 

PRINT  INPUT  DATA  B  79 

B  80 

WRITE  (G,117)  B  81 

WRITE  (B, 108)  (N, (C0RD(N,M).M=1,2),N=1,NP)  B  82 

WRITE  (G, 118)  B  83 

WRITE  (G,109)  (N, (N0P(N,M),M=1,4),IMAT(N),N=1»NE)  B  84 

WRITE  (B,119)  B  85 

WRITE  (6,110)  (NBC(I),NFIX(I),I=1,NB)  B  8B 

WRITE  (B, 120)  FDIN.FLD  B  87 

103  CONTINUE  B  88 

DO  104  IJ=1,200  B  89 

R10(IJ)=0.  B  90 

R20(IJ)=0.  B  91 

R30(IJ)=0.  B  92 

104  FORS(IJ)=0.  B  93 

DO  105  IJ=1,25  B  94 

105  NFIXK(IJ)=NFIX(IJ)  B  95 

RETURN  B  9B 

C  B  97 

lOG  FORMAT  (915)  B  98 

107  FORMAT  (13HONODAL  POINTS, 9X, I5/1X, 8HELEMENTS, 13X, IS/IX, 19HB0UNDARY  B  99 

1  CONDITIONS, 2X, I5/1X, 12HOUTPUT  LIMIT, lOX, I5/1X, 18HDEGREES  OF  FREED  B  100 

20M,3X,I5/1X,3HMATERIALS,12X,I5)  B  101 


108  FORMAT  (I10.2F10.3) 

109  FORMAT  (615) 

110  FORMAT  (215) 

111  FORMAT  (E10.3.2FIO.O) 

112  FORMAT  (I5.5F10.0) 

113  FORMAT  (3F10.0/I5.2F10.0) 

114  FORMAT  (215.2F10.0tI5) 

115  FORMAT  (I5.7X.3(F10.1.4X)<F5.3.7X.FB.S//') 

IIB  FORMAT  (1H1.I2.7A10) 

117  FORMAT  (14H0  MODAL  P0INTS/17X. 1HX» lOX. IHY) 

118  FORMAT  (lOHO  ELEMENTS/9X,1HI.4X.1HJ.4X.1HK»BX.3HMAT) 

119  FORMAT  (21H0  BOUNDARY  CONDITIONS) 

120  FORMAT  (IGHOPRINTING  SCHEME/5X.22H1.  REPORT  OUTPUT  EUERY.FB.2.2X.4 
IHMSEC/EX.REHR.  TERMINATE  OUTPUT  AT.FB.2.2X.4HMSEC) 

121  FORMAT  (1H0.20H  MATERIAL  PROPERTIES/IX.BHMAT.  N0..7X.2HE1. 12X.2HE2 
1, 11X.3HG12. 10X.3HU12, 10X.3HRH0/) 

122  FORMAT  (15H0BEAM  THICKNESS. 11X,F6.3/1X. lOHBEAM  WIDTH. 15X.FB.3^1X.  1 
14HSPHERE  DENSITY. 12X.F8.G/1X.  13HSPHERE  RADIUS. 12X.F6.3//1X. IIHIMPA 
2CT  NODE. 14X. I2/1X. 15HIMPACT  UELOCITY. lOX.FB.l/lX.  39HINTEGRATI0N  T 
3IME  INCREMENT!  X  E-OB  SEC).E10.3) 

123  FORMAT  (//.IX.  21HPERNANENT  DEFORMATION. 9X.F8.B) 

124  FORMAT  (/.IX.  15HUNL0ADING  POWER. 15X.FB.3) 


END 

SUBROUTINE  ESTIFM  (N) 

REAL  IB. LB 

COMMON  /CONTR/  TITLE(7).NPiNE.NB.NDr.NCN.NLD.NNAT.NSZF.LlrNT4.NDIN 
l.MATP.NPROB 

COMMON  /TIME/  T.DT.DDT.TAU.KCON.KCNT 
COMMON  /DISP/  01. 02.03.010.020.030 
COMMON  /DIMB/  TB.WB.PB.NQ.Dll 

COMMON  CORD(100.2).NOP(200.4).lNAT(2D0).ORT(25.S).NBC(25).NFIX(25) 
1.R1(200).R2(200).R3(200).R10(200).R2Q(200).R30(I200).FORS(200).SM(2 
200. 15) . SK(200. 15). ISP(200. 15) . SNPEM(200. 15) . ESTIF( 12. 12) . EMASS( 12. 
312).NFIXK(25) 

COMMON  /COMP/  0BR(3.3.25);ABD(S.6).TH(25).ZK(25)«MLAYER 
IB=WB»TB**3/12. 

LB=C0RD(N+1. l)-CORD(N.l) 

SOLB=LB*LB 

TPLB=LB*LB*LB 

1MN=IMAT(N) 

PARA1=0RT(IMN. l)*IB/70. 

IF  (MATP.EO.l)  PARAl=ABD(4.4)/70. 

ESTIF( 1 . 1 )=1200 . /TPLB*PARA1 
ESTIF( 1 .2)=B00  ./SQLB»PARA1' 

ESTtF( 1. 3)=30 . /LB«PARA1 
ESTIF(1.4)=-1200./TPLB*PARA1 
ESTIF  ( 1 . 5 )  =B00  .;/SOLB*PARA  1 
ESTIF  ( 1 . 6 )  =-30  ../LB*PARA1 
ESTIF(2.1)=ESTIF(1.2) 

ESTIF(2. 2)=384.!/LB*PARAl 
ESTIF(2.3)=22.»PARA1 
ESTIF( 2. 4 ) =-B0a . /S0LB*PARA1 
ESTIF(2,5)=21B.l/LB*PARAl  ' 

ESTIF(2,G)=-8.*PARA1 

ESTIF(3,1)=ESTIF(1.3) 

ESTIF(3.2)=ESTIF(2.3) 

ESTIF(3.3)=B.*LB»PARA1 

ESTIF(3.4)=-30.:/LB*PARAl 

ESTIF(3.5)=8.*PARA1 

ESTIF(3.B)=LB*PARA1 

ESTIF(4.1)=ESTIF(1.4) 

ESTIF(4.2)=ESTIF(2.4) 

ESTIF(4.3)=ESTIF(3.4) 

ESTIF(4.4)=120a./TPLB»PARAl 
ESTIF( 4 . 5 ) =-B0  0 . /SOLB*PAR A1 
ESTIF(4. B)*30 . /LB*PARA1 
ESTIF(5.1)=ESTIF(1.5) 

ESTIF(5.2)=ESTIF(2.5) 

ESTIF(5.3)=ESTIF(3.5) 


B  102 
B  103 
B  104 
B  105 
B  105 
B  107 
B  108 
B  109 
B  110 
B  111 
B  112 
B  113 
B  114 
B  115 
B  IIB 
B  117 
B  118 
B  119 
B  120 
B  121 
B  122 
B  123 
B  124 
B  125 
C  2 
C  3 
C  4 
C  5 
C  -B 
C  7 
C  8 
C  9 
C  10 
C  11 
C  12 
C  13 
C  14 
C  15 
C  16 
C  17 
C  18 
C  19 
C  20 
C  21 
C  22 
C  23 
C  24 
C  25 
C  26 
C  27 
C  28 
C  29 
C  30 
C  31 
C  32 
C  33 
C  34 
C  35 
C  36 
C  37 
C  38 
C  39 
C  40 
C  41 
C  42 
C  43 
C  44 
C  45 
C  46 
C  47 
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ESTIF(5»4)=ESTIF(4.5)  C  48 

ESTIF(5,5)=384./LB*PARA1  C  48 

ESTIF(5,6)=-S2.*PARA1  C  50 

ESTIF(G»1)=ESTIF(1»G)  C  51 

ESTIF(G,2)=ESTIF(2,G)  C  52 

ESTIF(G,3)=ESTIF(3,G)  C  53 

ESTIF(G,4)=ESTIF(4.G)  C  54 

ESTIF(G,5)=ESTIF(5,G)  C  55 

ESTIF(G»G)=G.»LB«PARA1  C  5G 

IF  (N.NE.l)  GO  TO  101  C  57 

WRITE  (G,103)  C  58 

WRITE  (G.102)  ((ESTIFCI.J),J=1»S)»I=1.6)  C  59 

101  CONTINUE  C  GO 

RETURN  C  G1 

C  C  G2 

102  FORMAT  (lX.GEll.B)  C  G3 

103  FORMAT  (1H0,38H  TYPICAL  STIFNESS  MATRIX  OF  AN  ELEMENT)  C  G4 

C  C  G5 

END  C  GG 

SUBROUTINE  EMASSM  (N)  D  2 

REAL  LB  D  3 

COMMON  /CONTR/  TITLE(7),NP,NE,NB»NDF.NCN,NLD»NMAT,NSZF,LI.NT4.NDIN  D  4 

l.MATP.NPROB  D  5 

COMMON  /TIME/  T, DT, DDT, TAU, KCON.KCNT  D  G 

COMMON  /DISP/  01,02,03,010,020,030  D  7 

COMMON  /DIME/  TB,WB,PB,N0,D11  D  8 

COMMON  CORD(100,2),NOP(200,4),IMAT(200),ORT(25,5),NBC(25),NFIX(25)  D  8 

1,R1(200),R2(200),R3(200),R10(200),R20(200),R30(200),FORS(200),SM(2  D  10 

200,15),SK(200,15),ISP(200,15),SMPEM(200,15),ESTIF(12,12),EMASS(12,  D  11 

312),NFIXK(25)  D  12 

COMMON  /COMP/  QBR(3,3,25),ABD(G,G),TH(25),ZK(25),MLAYER  D  13 

LB=C0RD(N+1,1)-C0RD(N, 1)  D  14 

AB=TB*WB  D  15 

SQLB=LB*LB  D  IG 

TPLB=LB»LB»LB  D  17 

QDLB=LB*LB*LB»LB  D  18 

IMN=IMAT(N)  D  19 

PARA2=ORT(IMN,5)»AB*LB/55440.  D  20 

EMASS(1,1)=21720.*PARA2  D  21 

EMASS(1,2)=3732.*LB*PARA2  D  22 

EMASS(1,3)=281.*SQLB*PARA2  D  23 

EMASS(1,4)=G000.*PARA2  D  24 

EMASS(1,5)=-1812.*LB*PARA2  D  25 

EMASS(1,G)=181.*S0LB*PARA2  D  2G 

EMASS(2,1)=EMASS(1,2)  D  27 

EMASS(2,2)=832.»SQLB»PARA2  D  28 

ENASS(2,3)=G3.*TPLB»PARA2  D  29 

EMASS(2,4)=1812.*LB»PARA2  D  30 

EMASS(2,5)=-532.*SQLB*PARA2  D  31 

EMASS(2,G)=52.*TPLB»PARA2  D  32 

EMASS(3,1)=EMASS(1,3)  D  33 

EMASS(3,2)=EMASS(2,3)  D  34 

EMASS(3,3)=G.*QDLB»PARA2  D  35 

EMASS(3,4)=181.*SQLB*PARA2  D  38 

EMASS(3,5)=-52.*TPLB*PARA2  D  37 

EMASS(3,G)=5o*QDLB*PARA2  D  38 

EMASS(4,1)=ENASS(1,4)  D  39 

EMASS(4,2)=EMASS(2,4)  D  40 

EMASS(4,3)=EMASS(3,4)  D  41 

EMASS(4,4)=21720.*PARA2  D  42 

EMASS(4,5)=-3732.*LB*PARA2  D  43 

EMASS(4,G)=281.*S0LB*PARA2  D  44 

EMASS(5, 1)=EMASS(1,5)  D  45 

EMASS(5,2)=EMASS(2,5)  D  4G 

EMASS(5,3)=EMASS(3,5)  D  47 

EMASS(5,4)=EMASS(4,5)  D  48 

EMASS(5,5)=832.*SQLB*PARA2  D  49 

EMASS(5,G)=-G9.*TPLB*PARA2  D  50 

EMASSCG, 1)=EMASS(1,G)  D  51 

EMASS(G,2)=EMASS(2,G)  D  52 


non  non  nnnnnnn  nnn  non  non  nnnn 
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EmSS(6,3)=EmSS(3,6)  D 

Ef1ASS(G,4)=EMASS(4,G)  D 

EMASS(G.5)=EmSS(5,G)  D 

EMASS(GfG)=G.*QDLB*PARA2  D 

IF  (N.riE.l)  GO  TO  101  D 

WRITE  (6,103)  D 

WRITE  (B, lOE)  ((EMASS(I,J), J=1,S),I=1,G)  D 

101  CONTINUE  D 

RETURN  D 

C  D 

102  FORMAT  (lX,6Eir.3)  D 

103  FORMAT  (1H0,34H  TYPICAL  MASS  MATRIX  OF  AN  ELEMENT)  D 

C  D 

END  D 

SUBROUTINE  FORMM  E 

E 

FORMS  MASS  MATRIX  E 

IN  UPPER  TRIANGULAR  FORM  E 

E 

COMMON  /CONTR/  TITLEC 7) , NP, NE, NB, NDF, NCN, NLD, NMAT, NSZF, LI, NT4, NDIN  E 
l,MATP,NPROB  E 

COMMON  /TIME/  T, DT,DDT,TAU,KCON, KCNT  E 

COMMON  /DISP/  01, 02, 03, 010,020,030  E 

COMMON  /DIMB/  TB, WB, PB, NO, Dll  E 

COMMON  CORD(100,E),NOP(E00,4),IMATC200),ORT(25,5),NBC(25),NFIX(E5)  E 
1,R1(2OO),R2(2OO),R3(2OO),R1O(20O),R2O(EOO),R3O(EOO),FORS(2OO),SM(2  E 
200,15),SK(200,15),ISP(200,15),SMPEM(200,15),ESTIF(12,1E),EMASS(12,  E 
312),NFIXKC25)  E 

COMMON  /COMP/  QBR(3,3,25),ABD(B,G),TH(25),ZK(25),MLAYER  E 

E 

SET  BANDMAX  AND  NO,  OF  EOUATIONS  E 

E 

NBAND=9  E 

E 

ZERO  MASS  MATRIX  E 

E 

DO  101  N=1,NSZF  E 

DO  101  M=1,NBAND  E 

101  SMCN,M)=0.  E 

E 

SCAN  ELEMENTS  E 

E 

DO  lOB  N=1,NE  E 

CALL  EMASSM  (N)  E 

E 

RETURNS  EMASS  AS  MASS  MATRIX  E 

E 

STORE  EMASS  IN  SM  E 

E 

FIRST  ROWS  E 

E 

DO  105  JJ=1,NCN  E 

NROWB=(NOP(N, JJ)-1)*NDF  E 

DO  105  J=1,NDF  E 

NR0WB=NR0WB+1  E 

I=(JJ-1)*NDF+J  E 

E 

THEN  COLUMNS  E 

E 

DO  104  KK=1,NCN  E 

NCOLB=  ( NOP  ( N,  KK )  )  *NDF  E 

DO  103  K=1,NDF  E 

L=(KK~1)*NDF+K  E 

NCOL=NCOLB+K+l“NROWB  E 

E 

SKIP  STORING  IF  BELOW  BAND  E 

E 

IF  (NCOL)  103,103,102  E 

102  SMC  NROWB , NCOL ) =SM ( NROWB, NCOL) +EMASS ( I ,  L )  E 

103  CONTINUE  E 


53 

54 

55 
5G 

57 

58 

59 
GO 
G1 
62 
G3 
B4 
65 
GG 

2 

3 

4 

5 
G 

7 

8 
9 

10 

11 

12 

13 

14 

15 
IG 

17 

18 

19 

20 
21 
22 

23 

24 

25 
2G 

27 

28 

29 

30 

31 

32 

33 

34 

35 
3S 

37 

38 

39 

40 

41 

42 

43 

44 

45 
4G 

47 

48 

49 

50 

51 

52 

53 

54 

55 
5G 
57 
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104  CONTIMUE 

105  CONTINUE 
lOG  CONTINUE 

C 

C  INSERT  BOUNDARY  CONDITIONS 

C 

DO  112  N=1.NB 
MX=10*«(NDF-1) 

I=NBC(N) 

NR0WB=(I-1)*NDF 

EXAMINE  EACH  DEGREE  OF  FREEDOM 

DO  111  M=1»NDF 
NR0WB=NR0WB+1 
ICON=NFIX(N)/NX 
IF  (ICON)  llOs 110» 10? 

SMCNROWB^ 1)=1. 

DO  lOS  J=2>NBAND 
SMCNROWB, J)=0. 

NR=NROWB+l-J 
IF  (NR)  109»109»108 
SM(NR. J)=0. 

CONTINUE 

NFIX(N)=NFIX(N)-NX*ICON 
NX=NX/10 
CONTINUE 

112  CONTINUE 

DO  115  N=1.NSZF 
K=0 

DO  114  M=1.NBAND 
MP=M-K 

IF  (ISP(N,M).LT.ISP(N»1))  GO  TO  113 
SM(N.MP)=SM(N,MP)+(DDT/G.)*SK(N,M) 
GO  TO  114 

113  K=K+1 

114  CONTINUE 

115  CONTINUE 

DO  IIB  I=1,NSZF 
DO  IIB  J=1.NBAND 
US  SMPEMCI. J)=SM(l9 J) 

c 

C  NRITE ( G>  1 )  ( ( Sfi  ( I  *  J ) » J=1  K  NBAMD ) » 1=1  •  NSZF ) 
C  1  FORMAT(2X.SE10.3) 

C 

RETURN 


C 

c 

c 


107 


108 

109 

no 

111 


c 


c 

c 

c 


c 

c 

c 


END 

COMMON^3raNTR/'^TITLE(7) .  NP,  NE,  NB,  NDF,  NCN.  NLD.  NMAT,  NSZF,  LI,.NT4,  NDIN 
L,MATP,NPROB 

COMMON  /TIME/  T,DT.DDT,TAU,KCON,KCNT 
COMMON  /DISP/  Ql, 02, 03, 010,020,030 

COMMON  CORD( 100?2)f NOP(200°4) , IMAT(200 ) , ORT (25, 5) , NBC(25) , NFIX(25) 

1,R1(200),R2(200),R3(200),R1D(200),R20(200),R30(200),FORS(200),SM(2 

200, 15),SK(200, 15), ISP(200, 15),SMPEM(200, 15),ESTIF(12, 12) ,EMASS( 12, 
312 ) 9 NF IXK ( 25 ) 

COMMON  /COMP/  QBR ( 3 » 3 >  25 ) ? ABD ( G » G ) >  TH( 25 ) >  ZK ( 25 ) » MLA YER 

SET  MAX.  NO.  OF  TERMS 

NMAX=0 

N0FF=9 


ZERO  ARRAYS 


DO  103  N=nNSZF 
DO  101  M=nNMAX 

SK(N9M)=0. 


E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

E 

F 

F 

F 

F 

F 

F' 

F 

F 

F 

F 

F 

F 

F 

F 

F 

F 

F 

F 

F 

F 

F 

F 


58 

59 
GO 
Bi 
G2 
G3 
G4 
G5 
6G 
67 
G8 
G9 

70 

71 

72 

73 

74 

75 
7B 

77 

78 

79 

80 
81 
82 

83 

84 

85 
8G 

87 

88 

89 

90 

91 

92 

93 

94 

95 
9G 

97 

98 

99 
100 
101 
102 

103 

104 

105 
2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 
IG 

17 

18 

19 

20 
21 
22 
23 


101 
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DO  102  M=2.N0FF  F  24 

102  ISP(N,M)=0  F  25 

103  ISP(N,1)=N  F  2G 

F  27 

SCAN  ELEMENTS  F  28 

F  29 

DO  no  N=1»NE  F  30 

CALL  ESTIFM  (N)  F  31 

F  32 

RETURNS  ESTIF  AS  STIFFNESS  MATRIX  F  33 

F  34 

STORE  ESTIF  IN  SK  WITH  A  TERM  IN  ISP  AS  A  POINTER  F  35 

F  3B 
F  37 

FIRST  THE  ROWS  F  38 

F  39 

1=0  F  40 

DO  109  JJ=nNCN  F  41 

NR0WB=(N0P(N.  JJ)-n*NDF  F  42 

DO  109  J=1.NDF  F  43 

NR0WB=NR0UB+1  F  44 

1=1+1  F  45 

F  48 

THEN  COLUMNS  OF  ESTIF  F  47 

F  48 

11=0  F  49 

DO  108  KK=1,NCN  F  50 

NC0LB=(N0P(N,KK)~1)*NDF  F  51 

DO  108  K=1»NDF  F  52 

NC0LB=NC0LB+1  F  53 

11=11+1  F  54 

F  55 

SEARCH  ISP  FOR  COLUMN  NO.  F  58 

F  57 

DO  105  M=1.N0FF  F  58 

IF  (ISP(NR0WB»M)-NC0LB)  104,107»104-  F  59 

104  IF  (ISP(NROWB.M))  108,108,105  F  60 

105  CONTINUE  F  81 

F  82 

FOUND  A  BLANK  NOW  STORE  NCOLB  F  83 

F  64 

108  ISP(NROWB,M)=NCOLB  F  65 

F  68 

NOW  STORE  ESTIF  F  87 

F  88 

107  SK(NR0WB,M)=ESTIF(1,II)+SK(NR0WB,M)  F  89 

F  70 

END  LOOP  ON  COLUMNS  F  71 

F  72 

108  CONTINUE  F  73 

F  74 

END  LOOP  ON  ROWS  F  75 

F  76 

109  CONTINUE  F  77 

F  78 

END  LOOP  ON  ELEMENTS  F  79 

F  80 

no  CONTINUE  F  81 

F  82 

INSERT  BOUNDARY  CONDITIONS  F  83 

F  84 

DO  114  N=1,NB  F  85 

NX=10**(NDF-1)  F  88 

I=NBC(N)  F  87 

NR0WB=(I-1)*NDF  F  88 

F  89 

EXAMINE  EACH  DEGREE  OF  FREEDOM  F  90 

F  91 

DO  113  M=1,NDF  F  92 

NROWB=NROWB+l  F  93 
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ICOM=NFIXK(N)/liX 
IF  (ICON) 

C 

C  STORE  ZERO  ON  DIAGONAL 

C 

111  SK(NROWB,1)=0.0 
NFIXK(N)=NFIXK(N)-‘NX*ICON 

112  NX=NX/10 

113  CONTINUE 

114  CONTINUE 
RETURN 


C 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


END 

cii  iRPni  ITTNF  LOAD 

comm  /CONTR/  title (?).  MP, me,  mb,  MDF,  MCN, MLB, MMAT,  MSZF,  LI ,  MT4, NDIM 
l,MATP,MPROB  __ 

COMMOM  /TIME/  T,  DT,DDT,TAU,KCOM,KCNT 
COMMON  /DISP/  Q1,Q2, 03, 010,020,030 
COMMOM  /DIMB/  TB, WB, PB, MO, D1 1 
COMMOM  /SPHERE/  STF,R,CABU(10),QKOM3T(10) 

COMMON  /PLASTIC/  DISPEM,MDISPEM,FORSPM,DISPM,QP 

COMMOM  CORD(100,2),NOP(200,4),IMAT(200),ORT(25,5),NBC(25),NFIX(25) 
1,R1(200),R2(200),R3(200),R10(200), R20 ( 200 ) , R30 ( 200 ) , f 2^^’ 

200,15),SK(200,15),ISP(200,15),SMPEM(200,15),ESTIF(12,12),EMASS(12, 

312)9  NF IXK ( 25 ) 

CONMON  /COMP/  QBR(393925)9 ABD(69G)»TH(25)9 ZK(25)9MLAYER 

SjPI=(4i/3i )isQRT(R)/( Cl .“'0RT(NMAT94)**2)/0RT(NMAT» l)+(l.-0RT(l94) 
^STFA=(4y/3! )*5QRT(R)/((1.-0RT(NMAT94)**2)/0RT(NMAT9 l)+l./0RT(l92)) 


STF=STFI 

IF  CMATP.EQ.l)  STF=STFA 
101  PAI=4.*ATAN(1. ) 

BALLM=(4./3. )#PAI*(R**3)*PB 


SIMPLY  SUPPORTED  SYMMETRY  CST1=0.5 
CLAMPED  CANTILEUER  CST1=1. 


CST1=1.0 


IFCNBCCl)  .EQ.  1)  CST1=0.5 


Q1=ACCEL.  OF  THE  MASS 
Q2=UEL0c  OF  THE  MASS 

q3=disp.  of  the  mass 


IF  CLI.GT.l.AND-KCNT.EQ.l)  GO  TO  102 
IF  (LI.GT.l.AND,KCNT.EQ-2)  GO  TO  103 
01=0.  . 


03=0. 


GO  TO  112 
102  010=01 
020=02 
030=03 

Q3=Q30+DT*Q20+0 . 5*DDT*Q10 
R3 ( NO ) =R30 ( NO ) +DT*R20 ( NO ) +0 
DIFD0=Q30“-R30(NQ) 
DIFDISP=03“R3(N0) 


o5*DDT*RlQ(NQ) 


WRITE(69400)  DIFDOtDIFDISP 

IF  (DIFDISP)  IlOf 1049 104 
1 03  Q3=Q30+DT*Q20+DDT*Q 1 0/3 . +DDT *01/6. 
DIFDO=030-R30(NO) 

DIFDISP=03“R3(NO) 


C 

C 

c 

c 


1RITE(G9400)  DIFDO, DIFDISP 

^00  F0RMATC/9  5X9?^DIFD0=?^9E15.39  5X9?JDIFDISP=?^tE15.3) 


F  94 
F  35 
F  96 
F  97 
F  98 
F  99 
F  100 
F  101 
F  102 
F  103 
F  104 
F  105 
F  106 
G  2 

G  3 

G  4 

G  5 

G  6 

G  7 

G  8 

G  9 

G  10 
G  11 
G  12 
G  13 
G  14 
G  15 
G  16 

G  17 

G  18 

G  19 

G  20 
G  21 
G  22 

G  23 

G  24 

G  25 

G  26 

G  27 

G  28 

G  29 

G  30 

G  31 

G  32 

G  33 

G  34 

G  35 

G  36 

G  37 

G  38 

G  39 

G  40 

G  41 

G  42 

G  43 

G  44 

G  45 

G  46 

G  47 

G  48 

G  49 

G  50 

G  51 

G  52 

G  53 

G  54 

G  55 

G  56 

G  57 

G  58 
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IF  (DIFDISP.lt. 0)  GO  TO  110  G  59 

104  IF  (DISPEM.EQ.0.0)  GO  TO  105  G  GO 

IF  ((DIFDISP.LT.DIFDO).AHD. (MDISPEM.EQ.O))  GO  TO  107  G  G1 

IF  ((DIFDISP.LT.DIFDO).AND. (NDISPEM.GT.O))  GO  TO  108  G  GS 

105  DO  lOG  J=1,NSZF  ,  G  G3 

lOG  FORS(J)=0.  '  G  G4 

F0RS(MQ)=STF*(DIFDISP)*»1.5»CST1  G  G5 

Q1=-F0RS(M0)/BALLM/CST1  G  GG 

IF  (KCMT.EQ.l)  GO  TO  113  G  G7 

IF  (KCNT.EQ.2)  GO  TO  109  G  G8 

107  NDISPEM=1  G  G9 

FORSPM=FORS(MQ)  G  70 

DISPM=DIFDO  G  71 

WRITE  (G, 114)  DISPEM.DISPM.DIFDISP.FORSPM  G  72 

IF  ((DIFDISP.LT.DISPEM).OR.(DISPM.LE.DISPEM))  GO  TO  111  G  73 

108  FORS(NQ)=FORSPM*((DIFDISP-DISPEM)/(DISPM-DISPEM))»»QP*CSTl  G  74 

Q1=-F0RS(NQ)/BALLM/CST1  G  75 

IF  (KCNT.EQ.l)  GO  TO  113  G  7G 

109  02=020+0. 5*DT*Q10+0.5*DT*Q1  G  77 

03=030+DT»020+DDT*010/3.+DDT»01/G.  G  78 

GO  TO  113  G  79 

no  FORS(liQ)=0.  G  80 

01=0.  G  81 

GO  TO  109  G  82 

111  LI=10000  G  83 

GO  TO  113  G  84 

112  FORS(NQ)=0.  G  85 

113  RETURN  G  8S 

C  G  87 

114  FORMAT  (///,5X,  7HDISPEM=, E10.3, 5X,  GHDISPM=.E10.3,5X.  8HDIFDIS  G  88 

1P=,E10.3,5X,  7HFORSPM=,E10.3)  G  89 

C  G  90 

END  G  91 

SUBROUTINE  HMTQ  H  2 

H  3 

SUBROUTINE  FOR  FINDING  (F)-(K)(U)  H  4 

H  5 

COMMON  /CONTR/  TITLE(7).NP»NE.NB,NDF.NCNf NLD.NMAT.NS2F»LI«NT4,NDIN  H  G 

l.MATP.NPROB  H  7 

COMMON  /TIME/  T,DT,DDT,TAU,KCON,KCNT  H  8 

COMMON  /DISP/  01,02.03,010,020.030  H  9 

COMMON  /DIMB/  TB.WB.PB.NO.Dll  H  10 

COMMON  CORD(100,2),NOP(200,4),IMAT(200),ORT(25,5),NBC(25),NFIX(25)  H  11 

1,R1(EOO),R2(200),R3(200),R10(200),R20(200),R30(200),FORS(200),SM(2  H  12 

200,15),SK(200, 15),ISP(200,15),SMPEM(200,15),ESTIF(12,12),EMASS(12,  H  13 

312),NFIXK(25)  H  14 

COMMON  /COMP/  0BR(3,3,25),ABD(B,G),TH(25),2K(25),MLAYER  H  15 

NT=9  H  IB 

DO  101  IJ=1,NSZF  H  17 

R1(IJ)=0.  H  18 

R2(IJ)=0.  H  19 

101  R3(IJ)=0.  H  20 

H  21 

H  22 

DO  105  N=1,NS2F  H  23 

FX=FORS(N)  H  24 

DO  102  M=1,NT  H  25 

L=ISP(N,M)  H  EG 

102  FX=FX-SK(N,M)*(R30(L)+DT»R2D(L)+(DDT/3.)*R10(L))  H  27 

IF  (SK(N,1))  104,103,104  H  28 

103  FX=0.  H  29 

104  R1(N)=FX  H  30 

105  CONTINUE  H  31 

RETURN  H  32 

H  33 

END  H  34 

SUBROUTINE  SOLUE  I  2 

I  3 

SPECIFICATION  STATEMENTS  I  4 

I  5 
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COMMON  /CONTR/  TITLE(7) , NP, NE, NB. NDF> NCN, NLD9 NMAT, NSZF9 LI ♦ NT49 NDIN 
1,MATP9NPR0B 

COMMON  /TIME/  T, DT, DDT, TAU, KCON, KCNT 
COMMON  /DISP/  01, Q2,Q3, 010,020,030 
COMMON  /DIME/  TB, WB, PB, NO, Dll 

COMMON  CORD(100,2),NOP(200,4), IMAT(200),ORT(25,5),NBC(25),NFIX(25) 
1,R1(200),R2(200),R3(200),R10(200),R20(200),R30(200),FORS(200),5M(2 
200,15),SK(200,15),ISP(200, 15) ,SMPEM(200, 15) , ESTIFC 12, 12), EMASSC 12, 
312),NFIXK(25) 

COMMON  /COMP/  QBRC3,3,25), ABD(B,G),TH(25),ZK(25),MLAYER 
NBAND:=9 

DO  101  I=i,NSZF 
DO  101  J=1,NBAND 

101  SMCI,  J)==SMPEM(I,  J) 

REDUCE  MATRIX 

DO  lOB  N=1,NSZF 
I=N 

DO  105  L=2,NBAND 
1=1+1 

IF  (SM(N,L))  102,105,102 

102  C=SM(N,L)/SM(N,  1) 

J=0 

DO  104  K=L,NBAND 
J=J+1 

IF  (SM(N,K))  103,104,103 

103  SM(I, J)=SM(I, J)“C^5M(N,K) 

104  CONTINUE 
SM(N,L)=C 

AND  LOAD  UECTOR 
FOR  EACH  EQUATION 

R1(I)=R1(I)-C*R1(N) 

105  CONTINUE 

lOG  R1(N)=R1(N)/SM(N, 1) 

BACK-SUBSTITUTION 

N=NSZF 

107  N=N-1 

IF  (N)  111,111,108 

108  L=N 

DO  no  K=2,NBAND 
L=L+1 

IF  (SM(N,K))  109,110,109 

109  R1(N)=R1CN)-SM(N,K)*R1CL) 
no  CONTINUE 

GO  TO  107 
111  RETURN 

END 

SUBROUTINE  INTEGTN 

COMMON  /CONTR/  TITLE (7) , NP, NE, NB, NDF, NCN, NLD, NMAT , NSZF , LI , NT4, NDIN 
l,MATP,NPROB 

COMMON  /TIME/  T, DT, DDT, TAU, KCON, KCNT 
COMMON  /DISP/  Q1,Q2,Q3,Q10,Q20,Q30 
COMMON  /DIME/  TB, WB, PB, NQ, Dll 

COMMON  CORD ( 100, 2) , NOP(200, 4 ) , IMAT(200) , DRTC25, 5) , NBCC25) , NFIX(25) 

1,R1(200),R2(200),R3C200),R10(200),REO(200),R30C200),FORS(200),SM(E 
200, 15),SK(200, 15), ISP(200, 15) , SMPEM(200, 15) , ESTIFC 12, 12) , EMASSC 12, 
312),NFIXKC25) 

COMMON  /COMP/  QBRC3,3,25),ABDCG,6),THC25),ZKC25),MLAYER 
COMMON  /PLOT/  NN, TTC25) , FF(25) , WC25) , UC25) 

R1=ACCEL.  OF  BEAM 
R2=UEL0,  OF  BEAM 
R3=DISPL.  OF  BEAM 
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17 

18 

19 

20 
21 
22 

23 

24 

25 
2G 

27 

28 

29 

30 

31 

32 

33 

34 

35 
3G 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 
43 

50 

51 

52 

53 

54 

55 
5G 

57 

58 
2 

3 

4 

5 
G 

7 

8 
9 

10 

11 

12 

13 

14 

15 
IG 

17 

18 


nnn  nnnn  non 


95 


DO  101  IJ=1,IHSZF  J  19 

R2(IJ)=R20(IJ)+0.5»DT«R10(IJ)+0.5*DT*R1(IJ)  J  20 

R3(IJ)=R30(IJ)+DT»R20(IJ)+(DDT/3.  )»R10(1J)  +  (DDT/G.  )*R1(IJ)  J  21 

101  CONTINUE  J  22 

IF  (KCNT.EQ.l)  GO  TO  107  J  23 

DO  102  1K=1,NP  J  24 

IK4=(IK-1)*3+1  J  25 

R30(IK)=R3(IK4)  J  2B 

102  CONTINUE  J  27 

IF  ((LI/10000). EQ.l)  GO  TO  103  J  28 

J  28 

PRINT  CONTROL  J  30 

J  31 

NTON=(LI-l)/NDIN  J  32 

IF  (NTON.NE.KCON)  GO  TO  105  J  33 

103  CONTINUE  J  34 

J  35 

SIMPLY  SUPPORTED  BEAM  CST2=2.  J  3B 

CANTILEUER  CST2=1.  J  37 

J  38 

CST2=1.  J  39 

J  40 

IF(NBC(1).EQ.1)CST2=2.  J  41 

J  42 

F=CST2*F0RS(NQ)  J  43 

APHA=Q3-R30(NQ)  J  44 

FF(NN)=F  J  45 

W(NN)=Q3*1000.  J  4G 

U(NN)=R30(NQ)*1000.  J  47 

WRITE  (G,108)  (TITLE(I), 1=1.7)  J  48 

T1=T*1.EG  J  48 

TT(NN)=T1  J  50 

NN=NN+1  J  51 

WRITE  (G,109)  T1.F.Q3,Q2.Q1.APHA  J  52 

DO  104  IK=1.NP  J  53 

IK3=IK*3  J  54 

STXX=R3(IK3)*TB/2.  J  55 

SIGX=ORT(l. 1)*STXX  J  5B 

STYY=-DRT(1.4)»STXX  J  57 

WRITE  (G.llO)  IK,R30(IK),STXX,STYY.SIGX  J  58 

104  CONTINUE  J  58 

KCON=KCON+l  J  GO 

105  CONTINUE  J  G1 

DO  lOB  IJ=1,NSZF  J  G2 

R10(IJ)=R1(IJ)  J  G3 

R20(IJ)=R2(IJ)  J  G4 

lOB  R30(IJ)=R3(IJ)  J  G5 

107  RETURN  J  SB 

C  J  G7 

108  FORMAT  (1H1.7A10///)  J  G8 

109  FORMAT  (10X.18HTIME  ELAPSED(MSEC). 13X.F7.3/10X.9HFORCE(LB).21X.E11  J  69 

1.3/10X,21HMASS  DISPLACEMENT(IN).9X,E11.3/10X,21HMASS  UELOCITY(IN/S  J  70 

2EC),9X,E11.3/10X,20HMASS  ACCEL. ( IN/SEC2). lOX. Ell .3/lOX, 15HINDENTAT  J  71 

3I0N(IN). 15X.E11.3///10X,4HNODE,9X,4HDISP. 13X. 9HSTRAIN-XX, 9X» 9HSTRA  J  72 

4IN-YY.9X.3HSTRESS-XX/)  J  73 

110  FORMAT  (9X.I3.7X,E12.3,7X.E12.3,7X.E12.3.7X.E12.3)  J  74 

C  J  75 

END  J  7G 

SUBROUTINE  CMPD  K  2 

COMMON  /CONTR/  TITLE(7),NP,NE,NB.NDF.NCN,NLD,NMAT,NSZF,LI.NT4.NDIN  K  3 

l.MATP.NPROB  K  4 

COMMON  /TIME/  T. DT, DDT. TAU. KCON, KCNT  K  5 

COMMON  /DISP/  01,02.03.010.020.030  KG 

COMMON  /DIMB/  TB.WB.PB.NQ.Dll  K  7 

COMMON  CORD(100.2),NOP(200,4),IMAT(200),ORT(25.5),NBC(25).NFIX(25)  K  8 

1,R1(200),R2(200),R3(200),R10(200),R20(200),R30(200).FORS(200),SM(2  K  9 

200,15),SK(200. 15), ISP(200,15),SMPEM(200.15),ESTIF(12.12).EMASS(12.  K  10 

312),NFIXK(25)  K  11 

COMMON  /COMP/  QBR(3.3.25),ABD(G,B).TH(25).ZK(25),MLAYER  K  12 

DIMENSION  0(3,3),  TK(25)  K  13 
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DO  102  J=1.3  K  14 

DO  102  K=l,3  K  15 

ABD(J+3»K+3)=0.0  K  IG 

ABD(J+3.K)=flBD(J+3,K+3)  K  17 

ABD(J,K+3)=ABD(J+3,K)  K  18 

ABD(J»K)=ABD(J,K+3)  K  19 

DO  101  1=1,25  K  HO 

QBRCJ,K, I)=0.  K  21 

101  CONTIMUE  K  22 

102  CONTINUE  K  23 

READ  (5,108)  MLAYER  K  24 

M=MLAYER  K  25 

READ  (5,109)  (L,THCL),TK(L), 1=1, M)  K  2B 

TTK=0.0  K  27 

ZK(1)=TTK  K  28 

DO  103  1=1, M  K  28 

TTK=TTK+TK(I)  K  30 

ZK(I+1)=TK(I)+ZK(I)  K  31 

103  CONTINUE  l<  32 

MM=M+1  K  33 

DO  104  1=1, MM  K  34 

ZK(I)=ZI<(I)-TTK/2.  K  35 

104  CONTINUE  K  3B 

DEL=4.*ATAN(1. )/180.  K  37 

DEN=l.-0RT(l,2)*0RT(l,4)«-»2/0RT(l,l)  K  38 

Q(l,l)=ORT(l,l)/DEN  K  39 

Q(2,2)=0RT(1,2)/DEN  K  40 

Q(2,1)=0RT(1,4)»Q(2,2)  K  41 

Q(1,2)=Q(2,1)  K  4H 

Q(3,3)=0RT(1,3)  K  43 

Q(3,2)=0.0  K  44 

Q(3,1)=Q(3,2)  K  45 

Q(2,3)=Q(3,1)  K  4B 

Q(1,3)=Q(2,3)  K  47 

DO  105  1=1, M  K  48 

ANGL=TH(I)»DEL  K  49 

C=COS(ANGL)  K  50 

W=SIN(ANGL)  K  51 

QBR(1,1,I)=Q(1,1)-»C**4+2.*(Q(1.2)+2.*Q(3»3))»(C*W)**2+Q(2,2)*W»  K  52 

1  «4  K  53 

QBR(2,1, I)=(Q(1, 1)+Q(2,2)-4.*Q(3,3))*(C*M)**2+Q(1,2)»(W**4+C**4  K  54 

1  )  K  55 

QBR(1,2,I)=QBR(2,1,I)  K  5B 

QBR(2,2,I)=Q(1,1)»U*»4+2.*(Q(1,2)+2.*Q(3,3))»(C»U)«*2+Q(2,2)*C*  K  57 

1  *4  K  58 

QBR(3,1, 1)  =  (Q(1,  l)-Q(l,E)-2.*Q(3,3))*W*C«*3+(Q(l,H)-Q(2,2)+(2.)  K  59 

1  *0(3,3) )«(W)*(C**3)  K  60 

QBR(1,3,I)=QBR(3,1,I)  K  G1 

QBR(3,2, I)=(Q(1, 1)-Q(1,2)-2.*Q(3,3))»W**3*C+(Q(1,2)-Q(2,2)+2.*Q  K  SB 

1  (3,3))*U«C«*3  K  G3 

QBR(2,3,I)=QBR(3,2,I)  K  B4 

QBR(3,3, I)=(Q(1, 1)+Q(B,B)-2.*Q(1,2)-2.*Q(3,3))*(W«C)*»2+Q(3,3)*  K  G5 

1  (U**4+C**4)  K  BB 

K  B7 

WRITE(B,500)  I,TH(I),TI<(I)  K  G8 

WRITE(B,510)  K  G9 

WRITE(G,5E0)  ( (QBR(J,K, I),K=1,3), J=l,3)  K  70 

K  71 

105  CONTINUE  K  72 

DO  107  J=l,3  K  73 

DO  107  K=l,3  K  74 

DO  lOG  1=1, M  K  75 

ABD(J,K)=ABD(J,K)+QBR(J,K,I)*(ZK(I+1)-ZK(I))  K  7G 

ABD(J,K+3)=ABD(J+3,K)+QBR(J,K,I)*(ZK(I+l)**2-ZK(I)*»2)/2.  K  77 

ABD(J+3,K)=ABD(J,K+3)  K  78 

ABD ( J+3 , K+3 ) = ABD ( J+3 ,  l<+3 ) +QBR ( J , K , I ) * ( ZK ( 1+ 1 ) «*3-ZK ( I ) **3 ) /3  K  79 

1  .  K  80 

lOB  CONTINUE  K  81 

107  CONTINUE  K  82 

WRITE  (G,110)  K  83 
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WRITE  (G.lli)  C(ABD(I»J),J=1.G),I=1»G)  K  84 

K  85 

500  F0RMAT(2X,«LAYER=»,ia,5X,»ANGLE=*,F5.2,5X.*THICKNESS=»»F7.3)  K  8B 

510  F0RMAT(2X.*QBAR-I1ATRIX*)  K  87 

520  F0RMAT(5X.3E12.3/)  K  88 

K  80 

RETURN  K  30 

C  K  91 

108  FORMAT  (15)  K  02 

109  FORMAT  (I5.F5.0.F10.0)  K  93 

110  FORMAT  (///.lX.lOHABD  MATRIX)  K  94 

111  FORMAT  (1X.GE11.3)  K  95 

C  K  9G 

END  K  97 
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APPENDIX  B 

A  COMPUTER  PROGRAM  FOR  ESTIMATING  THE  CONTACT  FORCE  HISTORY  BY  USING  THE 
EQUIVALENT  MASS  MODEL 

This  program  has  been  written  for  simply-supported  beams  only. 
Cantilever  beams  and  simply-supported  plates  will  be  added  in  the 
near  future. 

This  program  will  be  a  subprogram  in  a  large  finite  element 
program  capable  of  analyzing  impact  responses  of  beams  and  plates. 

This  subprogram  will  be  used  to  provide  an  estimate  of  the  contact 
time  so  that  one  may  select  a  proper  time  increment  for  the  finite 
difference  used  in  the  program.  For  this  reason,  the  input  cards 
for  this  subprogram  were  written  to  be  identical  to  that  for  the 
program  presented  in  Appendix  A. 
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PROGRAM  MAIN  (INPUT, OUTPUT, TAPE5=INPUT,TAPEG=0UTPUT)  A  2 

COMMON  /CONTR/  TITLE( 10 ) , NP, NE, NB, NDF, NCN, NLD, NMAT, NSZF, LI, NT4, ND+  A  3 

1N,MATP,NPR0B  A  4 

COMMON  /TIME/  T, DT, DDT, TAU, KCON, KCNT  A  5 

COMMON  /DISP/  Q1,Q2, 03, 010,020,030  A  G 

COMMON  /DIME/  TB, WB, PB, NO, Dll  A  ? 

COMMON  /SPHERE/  STF, R, CABU(10),QKONST( 10)  A  8 

COMMON  CORDC100,2),NOP(200,4),IMAT(200),ORT(25,5),NBC(25),NFIXC25)  A  9 

1,R1(200),R2(200),R3(200),R10(200),R2D(200),R30(200),FORS(200),SM(2  A  10 

200,15),SK(200,15),ISP(200,15),SMPEM(200,15),ESTIF(12,12),EMASS(12,  A  11 

312),NFIXK(25)  A  12 

COMMON  /COMP/  QBR(3,3,25),ABD(6,G),TH(25),ZK(25),MLAYER  A  13 

COMMON  X1,X2,ND1,ND2  A  14 

REAL  IB  A  15 

A  IG 

READ  AND  PRINT  TITLE  AND  CONTROL  A  17 

A  18 

101  READ  (5,117)  NPROB, (TITLE(I),I=1,10)  A  19 

IF  (NPROB. EQ.O)  GO  TO  107  A  20 

WRITE  (G,118)  NPROB, (TITLE(I), 1=1, 10)  A  21 

READ  (5,108)  NP,NE,NB,NTM,NMAT,NDIN,MATP,NDC,I1  A  22 

NDF=3  A  23 

NLD=NDIN*NTM  A  24 

FLD=FLOAT(NLD)/10.  A  25 

FDIN=FLOAT(NDIN)/10.  A  2G 

READ  (5,114)  TB,WB,R,NQ,QE,DT  A  27 

WRITE  (G, 109)  NP,NE,NB, NLD, NDF, NMAT  A  28 

NLD=NLD+1  A  29 

A  30 

READ  AND  PRINT  MATERIAL  DATA  A  31 

SPHERE  DATA!  L=NMAT  (LAST  NAT.  CARD)  A  32 

A  33 

READ  (5,113)  (L, (ORT(L, I), 1=1,5), N=l, NMAT)  A  34 

PB=0RT(NMAT,5)  A  35 

WRITE  (6,123)  TB,WB,PB,R,NQ,Q2,DT  A  3G 

WRITE  (6,122)  A  37 

WRITE  (G, 116)  (N, (ORT(N, I), 1=1,5), N=l, NMAT)  A  38 

A  39 

READ  INDENTATION  CARD  A  40 

A  41 

READ  (5,112)  STF  A  42 

A  43 

READ  NODAL  POINT  DATA  A  44 

AND  A  45 

READ  ELEMENT  DATA  A  4G 

A  47 

DO  103  1=1, NDC  A  48 

READ  (5,115)  ND1,ND2,X1,X2,IMT  A  49 

EL=(X2-X1)/FL0AT(ND2~ND1)  A  50 

C0RD(ND1,1)=X1  A  51 

C0RD(ND2,1)=X2  A  52 

CORD(ND2,2)=0.0  A  53 

C0RD(ND1,2)=C0RD(ND2,2)  A  54 

NDD=ND2-1  A  55 

DO  102  J=ND1,NDD  A  5G 

CORD(J+l,l)=CORD(J,l)+EL  A  57 

CORD(J+1,2)=0.0  A  58 

NOP(J,l)=J  A  59 

N0P(J,2)=J+1  A  GO 

NOP(J,4)=0  A  61 

N0P(J,3)=N0P(J,4)  A  62 

IMAT(J)=IMT  A  63 

102  CONTINUE  A  64 

103  CONTINUE  A  65 

A  66 

READ  BOUNDARY  DATA  A  67 

A  68 

READ  (5,111)  (NBC(I),NFIX(I),I=1,NB)  A  69 

IF  (MATP.EQ.l)  CALL  CMPD  A  70 

C  A  71 
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ISOTROPIC 

COMPOSITE 


MATP=0.0 
MATP=1 . 0 


IF  (Il.NE.O)  GO  TO  104  " 

PRINT  INPUT  DATA  ^ 

P) 

WRITE  (G.llO)  (Nj (C0RD(N»M).M=1»2).N=1»NP)  ^ 

A 

WRITECg! 3) (N» (N0P(N. M) , M=1. 4) . IMATCN) , N=l. NE)  ^ 

WRITE  (G,120)  A 

WRITE  (G>111)  (NBC(I)>NFIX(I)» I=1»N3) 

WRITE  (G,121)  FDIN.FLD  2 

104  CONTINUE  p, 

DO  105  IJ=1>200  A 

R10(IJ)=0.  A 

R20(IJ)=0.  A 

R30(IJ)=0.  A 

105  F0RS(IJ)=0.  A 

DO  lOG  IJ=1.25  A 

lOG  NFIXK(IJ)=NFIX(IJ)  ^ 

CALL  TMX  P) 

3  F0RI1AT(GI5)  ^ 

103  FORMATdOHO  ELEMENTS/9X»  1HI*4X*  1HJ»4X»  1HK»8X»3HMAT)  A 

A 

GO  TO  101  PI 

107  STOP  p, 

PI 

ini  M^HnNOIlAL  POINTS,  SX.  I5/lX.8HELEt1ENTS.  13X,  I5/1X.  19HB0UNDARY  A 

G|1IT,10X,I5^1X,18HDEGREES  OF  FREED  0 
S0M,3X,I5/1X,9HMATERIALS,12X,I5)  " 

110  FORMAT  (I10,2F10.3)  P, 

111  FORMAT  (215)  p, 

112  FORMAT  (E10.3)  p, 

113  FORMAT  (I5,5F10.0)  p, 

114  FORMAT  (3F10.0/I5,2F10.0) 

115  FORMAT  (2I5,2F10.0,I5)  p, 

IIG  FORMAT  (I5,7X,3(F10.1,4X),F5.3,7X,F8.B//) 

117  FORMAT  (I2,10A7)  p, 

119  format  (HHO^NODAL^POINTS/ITX.IHX.IOX.IHY)  ^ 

ii?  =  EUERV,FG.2.2X,4  A 

122^Sf  aH??i0H“A[  2 

l?3^™MAT^(liH0BEAM^THICKNEsi(IN^^^  14HBEAM  WIDTH(IN),  11X,FG  A 

123  FORMAT  (13H0BEAMTHICKNE5bU^  l^HSPHERE  RADIUSCIN) , 8X  A 

LIg  3//?x?UHIMpScrSDErU?I2/lX.^^^  UELOCITYC IN/SEC), 2X.F  A 

Ig.PiX^  39HINTEGRATI0N  TIME  INCREMENTC  X  E-OG  SEC).E10.3)  A 

:  A 

END  B 

COMMON^/CONTR/  T I TLE ( 1 0 ) ,  NP , NE » NB , NDF , NCN , NLD , NMAT , NSZF , LI , NT4 , ND I  B 

1N,MATP,NPR0B  r 

COMMON  /TIME/  T,DT,DDT,TAU,KCON,KCNT  ° 

COMMON  /DISP/  01,02,03,010,020,030  ^ 

Sf!Suo!.“S'(™i?i!.imT<EOO),ORT<E5.5).NBC.aS>,NFIX<E5>  B 

1  PI  I’pnni  p?fPf)0),R3(200),R10(200),R20(200),R30(200),FORS(200),SM(2  B 

E05di??sS!5)dSP?20O,l5)7sSpEM(EOO,15),ESTIF(  B 

^COMMON  /COMP/  QBR(3,3,25),ABD(G,G),TH(25),ZK(E5),MLAYER  B 

COMMON  X1,X2,ND1,ND2 


nnnnn 


101 


REAL  IB 

DIMENSION  Q(3,3),  TK(25) 

DO  102  J=l,3 
DO  102  K=l,3 

ABD(J+3fK+3)=0.0 

ABD ( J+3 , K ) = ABD ( J+3 • K+3 ) 

ABD(J,K+3)=ABD(J+3,K) 

ABD(J,K)=ABD(JfK+3) 

DO  101  1=1,25 
QBRCJ,K,I)=0. 

101  COMTIMUE 

102  CONTINUE 

READ  (5,108)  MLAVER 
M=MLAYER 

READ  (5,109)  (L,TH(L),TK(L),I=1,M) 

TTK=0,0 
ZK(1)=TTK 
DO  103  1=1, M 
TTK=TTK+TK(I) 

ZK(I+l)=TK(r)+ZK(I) 

103  CONTINUE 
MM=M+1 

DO  104  1=1, MM 

ZK(I)=ZK(I)“TTK/2. 

104  CONTINUE 
DEL=4.*ATAN(1.)/180. 

DEN= 1 . -ORT ( 1 , 2 ) *ORT (1,4) **2/0RT (1,1) 

Q(l,l)=ORT(l,l)/DEN 

Q(2,2)=0RT(1,2)/DEN 

Q(2,1)=0RT(1,4)*Q(E,2) 

Q(1,2)=Q(2,1) 

Q(3,3)=0RT(1,3) 

Q(3,2)=0.0 

Q(3,1)=Q(3,2) 

Q(2,3)=Q(3,1) 

Q(1,3)=Q(2,3) 

DO  105  1=1, M 

ANGL=TH(I)*DEL 

C=COS(ANGL) 

W=SIN(ANGL) 

QBR(1,1,I)=Q(1,1)*C*»4+2**(Q(1,2)+2.*Q(3,3))*(C*W)»*2+Q(2,2)*W» 
1  *4 

QBR(2,1,I)=(Q(1,1)+Q(2,2)-4.^^Q(3,3))»(C*U)«*2+Q(1,2)»(W**4+C**4 
1  ) 

QBR(1,2,I)=QBR(2,1,I) 

QBR(2,2, I)=Q(l,l)*W**4+2,*(Q(l,2)+2.*Q(3,3))*(C<tW)*»2+Q(2,2)*C* 
1  *4 

QBR(3,l,I)  =  (Q(l,l)-Q(l,2)“2.*Q(3,3))*W‘:^C**3+(Q(l,2)-'Q(2,2)  +  (2.) 
1  *Q(3,3))*(W).*(C**3) 

QBR(1,3,I)=QBR(3,1,I) 

QBR(3,2,I)=(Q(1,1)“Q(1,2)-2.*Q(3,3))*W**3*C+(Q(1,2)“Q(2,2)+2.*Q 
1  (3,3))*W*C**3 

QBR(2,3,I)=QBR(3,2,I) 

QBR(3,3, I)=(Q(1,1)+Q(2,2)“2.»Q(1,2)“2.*Q(3,3))*(W*C)**2+Q(3,3)* 

1  (|^**4+Q**4) 

WRITE(G,500)  I,TH(I),TK(I) 

WRITE(G,510) 

URITE(G,520)  ( (QBR( J, K, I ) , K=l, 3) , J=l, 3) 

105  CONTINUE 

DO  10?  J=l,3 
DO  10?  K=l,3 
DO  lOG  1=1, M 

ABD(J,K)=ABD(J,K)+QBR(J,K,I)*(ZK(I+1)“ZK(.I)) 

ABD(J,K+3)=ABD(J+3,K)+QBR(J,K,I)*(ZK(I+l)**2-ZK(I)**2)/2.. 

ABD(J+3,K)=ABD(J,K+3) 

ABD ( J+3 , K+3 ) = ABD ( J+3 , K+3 ) +QBR ( J , K , I ) * ( ZK ( I + 1 ) **3-ZK ( I ) **3 ) /3 

1 

106  CONTINUE 


B  14 
B  15 
B  IG 
B  1? 
B  18 
B  19 
B  20 
B  21 
B  22 
B  23 
B  24 
B  25 
B  2G 
B  2? 
B  28 
B  29 
B  30 
B  31 
B  32 
B  33 
B  34 
B  35 
B  3G 
B  3? 
B  38 
B  39 
B  40 
B  41 
B  42 
B  43 
B  44 
B  45 
B  4G 
B  4? 
B  48 
B  49 
B  50 
B  51 
B  52 
B  53 
B  54 
B  55 
B  5G 
B  5? 
B  58 
B  59 
B  GO 
B  G1 
B  62 
B  G3 
B  G4 
B  G5 
6  GG 
B  ^  G? 
B  G8 
B  69 
B  ?0 
B  ?1 
B  ?2 
B  ?3 
B  ?4 
B  ?5 
B  ?G 
B  ?? 
B  ?8 
B  ?9 
B  80 
B  81 
B  82 
B  83 
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102 


107  CONTINUE 

WRITE  (G,lll)  ((ABDCI, J).J=1»6),I=1,B) 

WRITE  (G,110) 

500  FORMAT ( 2>U  *LAYER=^f , IH. 5X, *ANGLE=*»  F5. 2, 5Xf  *THICKNESS=*f  F7. 3 ) 
510  FORMAT (2X,*QBAR-MATRIX*) 

520  F0RMAT(5Xf3E12*3/) 

RETURN 


108  FORMAT  CI5) 

109  FORMAT  (I5,F5,0.F10.0) 

no  FORMAT  (///nXnOHABD  MATRIX) 
111  FORMAT  (lX,6En.3) 


END 

SUBROUTINE  TMX  ^ 

COMMON  /CONTR/  T ITLE ( 1 0 ) >  NP, NE»  NBr  NDF^  NCN >  NLDf  NMAT » NSZF t  LI , NT4 , NDI 
1N,MATP,NPR0B 

COMMON  /TIME/  T. DTt DDT, TAU, KCON, KCNT 
COMMON  /DISP/  01,02,03, 010,020,030 
COMMON  /DIMB/  TB, WB, PB, NO, D1 1 
COMMON  /SPHERE/  STF, R, CABUC 10) ,OKONST(10) 

COMMON  /PLASTIC/  DISPEM, NDISPEM, FORSPM, DISPM 

COMMON  CORDC 100, 2), NOP (200,4 ), IMAT(200) ,ORT (25,5), NBC(25),NFIX(25) 
1,R1C200),R2(200),R3(200),R10(200),R20(200),R30(200),FORS(200),5M(2 
200,15),SK(200,15),ISP(200,15),SMPEMC200, 15) , ESTIFC 12, 12) , EMASSC 12, 
312),NFIXK(25) 

COMMON  /COMP/  QBR(3,3,25),ABD(6,G),TH(25),ZK(25),MLAYER 
COMMON  X1,XE,ND1,ND2 
REAL  IB 


EQUIUALENT  MODEL  FOR  ESTIMATING  TIME  INTERUAL 


IF  (STF. NE. 0.0)  GO  TO  101 

STFI=(4./3. )*SQRT(R)/((1.-0RT(NMAT,4)**2)/0RT(NMAT, 1)+(1.~0RT(1,4) 
1**2)/0RT(1,  D) 

STFA=(4./3. )*SQRT(R)/((1.-0RT(NMAT,4)**2)/0RT(NMAT, 1)+1./0RT(1,2)) 


STF=STFI 

IF  (MATP.EQ.l)  STF=STFA 
101  PAI=4o*ATAN(l.) 

BALLM= ( 4 . /3 . ) *P A I * ( R**3 ) *PB 


BL=X2”X1 

AB=WB*TB 

IB=WB*TB*^3/12. 

WATP=FLOAT(MATP) 

Dll=ORT(i,l)^IB 
IF  (MATP.EQ.l)  D11=ABD(4,4) 

WRITE  (6,105)  BL,STF,WATP,D11,ABD(4,4) 
WN1=( (D11*(PAI/BL)**4)/(AB^0RT(1,5) ) )**0.5 


TN1=2.*PAI/WN1 

TD=0.0 

DO  103  1=1, 1000 


N=0 

TD=TD+n0E“-7 


F1=0.0 


G1=0.0 

DO  102  J=l,100 


N=N+1 

C=FL0AT ( NQ“ 1 ) /FLOAT ( ND2~ND1 ) 
WX=SIN(N*PAI*C) 

PP=1 .  /  ( 4 ,  *N*^^4*TD**2-TN1**2 ) 

QQ=TN 1/(2. *N»*2*TD ) 

SS=WNl*N*^2*TD/2. 

SUMF=(PP*N*^^2*(  1~QQ*SIN(SS)  )»WX)«*2 
SUMG= ( PP^COS ( SS ) *WX ) **2 
F1=F1+5UMF 
G1=G1+SUMG 
CONTINUE 

SUMl=2.*Fl*16.*BL*^f3*TD**2/(Dll*PAI**2) 
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27 

28 
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39 

40 

41 
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103 


SUM2=2.*Gl*lG.*0RT(l,5)*(=iB*BL**7/(Dll**2*PftI**4)  C  56 

EMT=1./(SUM1+SUM2)  C  57 

STFE=1./EMT+1./BALLM  C  58 

STFT=STF*STFE  C  59 

APHAMAX=(1.25*Q2**2/STFT)**0.4  C  60 

T0TALT=2.94*APHAMAX/Q2  C  61 

FriAX=STF*APHAMAX**1.5  C  62 

EPS=TOTALT-TD  C  63 

IF  (ABS(EPS).LE.l.E-7)  GO  TO  104  C  64 

103  COMTINUE  C  65 

104  WRITE  (6,106)  TOTALT,TD,F(1AX, APHAMAX,EI1T  C  66 

RETURN  C  67 

C  C  G8 

105  FORMAT  (//,5X,  5HBEAM=,E10.3,5X,  4HSTF=,E10.3,5X,  5HMATP=,E10.3  C  69 

1,5X.  4HD11=,E10.3,5X,  9HABD(4, 4)=, E10.3)  C  70 

106  FORMAT  (//,5X,  7HTOTALT=, E10.3, 5X,  3HTD=,E10.3,5X,  5HFMAX=.E10.  C  71 

13,5X,  8HAPHAMAX=,E10.3,5X,  4HEMT=, E10.3)  C  72 

C  C  73 

END  C  74 
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